Clase 5 Más sobre problemas de clasificación

En problemas de clasificación, queremos usar información del modelo para tomar cierta acción. Ejemplos típicos de que requieren de resolver un problema de clasificación son:

• Dar un tratamiento a una persona para combatir una enfermedad. El tratamiento tienen costos monetarios y efectos secundarios. ¿Cuándo deberíamos tratar o hacer seguimiento de una persona?
• Decidir si hacemos un descuento a un cliente que tiene probabilidad alta de cancelar su contrato en los siguientes 3 meses.
• Para una búsqueda dada de restaurantes, decidir qué restaurantes debemos poner en las primeras posiciones de los resultados de la búsqueda.
• Decidir si una imagen contiene una persona o no, con el fin de activar una alarma en ciertas condiciones.

En la mayoría de estos ejemplos, no queremos encontrar un clasificador si un cliente va abandonar o no, tiene una enferdad o no o si un segmento de imagen contiene una persona o no, etc. Esto solo aplica para problemas con tasas de ruido muy bajas donde podemos separar claramente las clases - lo cual no es tan común, especialmente en problemas de negocios, por ejemplo.

Igual que en regresión producir intervalos de predicción (en problemas que no son de ruido bajo) permiten tomar mejores decisiones downstream del modelo, en clasificación producir probabilidades de clase permite tomar mejores decisiones que toman en cuenta aspectos del problema particular que nos interesa.

En todos los problemas de arriba, la dificultad es que al tomar la decisión de clasificar un caso un una clase específica, para los cuales se va a llevar a cabo una acción, diversos costos intervienen cuando cometemos distintos errores:

• Por ejemplo, diagnosticar a alguien con una enfermedad cuando no la tiene tiene consecuencias distintas a diagnosticar como libre de enfermedad a alguien que la tiene. Estas consecuencias dependen de cómo son son los tratamientos consecuentes, y de qué tan peligrosa es la enfermedad.

• Cuando usamos un buscador como Google, es cualitativamente diferente que el buscador omita resultados relevantes a que nos presente resultados irrelevantes.

En general, los costos de los distintos errores son distintos, y en muchos problemas quiséramos entenderlos y controlarlos individualmente. Aunque en teoría podríamos asignar costos a los errores y definir una función de pérdida apropiada, en la práctica esto muchas veces no es tan fácil o deseable.

Cuando producimos salidas que son clasificadores “duros” (asignan a una clase, por ejemplo, usando máxima probabilidad u otro método) varios problemas pueden aparecer. El desempeño de clasificadores duros generalmente se mide con variación de la matriz de confusión:

Matriz de confusión. Sea \(\hat{G}\) un clasificador. La matriz de confusión \(C\) de \(\hat{G}\) está dada por \(C_{i,j} =\) número de casos de la clase verdadera \(j\) que son clasificados como clase \(i\) por el clasificador

Ejemplo

En un ejemplo de tres clases, podríamos obtener la matriz de confusión:

	A	B	C
A.pred	50	2	0
B.pred	20	105	10
C.pred	20	10	30

Esto quiere decir que de 90 casos de clase \(A\), sólo clasificamos a 50 en la clase correcta, de 117 casos de clase \(B\), acertamos en 105, etcétera. Podemos ver esta tabla de distintas formas, por ejemplo, usando porcentajes por columna, nos dice cómo se distribuyen los casos de cada clase:

knitr::kable(round(prop.table(tabla_1, 2), 2))

	A	B	C
A.pred	0.56	0.02	0.00
B.pred	0.22	0.90	0.25
C.pred	0.22	0.09	0.75

Mientras que una tabla de porcentajes por renglón nos muestra qué pasa cada vez que hacemos una predicción dada:

knitr::kable(round(prop.table(tabla_1, 1), 2))

	A	B	C
A.pred	0.96	0.04	0.00
B.pred	0.15	0.78	0.07
C.pred	0.33	0.17	0.50

Ahora pensemos cómo podría sernos de utilidad esta tabla. Discute:

• El clasificador fuera uno de severidad de emergencias en un hospital, donde A=requiere atención inmediata B=urgente C=puede posponerse (poco adecuado).

• El clasificador fuera de tipos de cliente de un negocio. Por ejemplo, A = cliente de gasto alto, B=cliente medio, C=cliente de gasto bajo. Tenemos un plan para incrementar la satisfacción de los clientes: para clientes de gasto bajo cuesta muy poco, para los clientes de gasto medio tiene precio bajo, y cuesta mucho para los clientes de gasto alto (más adecuado).

La tasa de incorrectos es la misma en los dos ejemplos, pero la adecuación del clasificador es muy diferente.

Nótese que un clasificador bueno, en general, es uno que tiene la mayor parte de los casos en la diagonal de la matriz de confusión. Es difícil decir entonces cuándo un clasificador “duro” es bueno o malo sin tener más datos acerca del problema (a menos que su tasa de incorrectos sea cercana a 0).

5.1 Ejemplo: decisiones basadas en probabilidades

Supongamos que tenemos un plan para retener a clientes. Construimos un modelo que nos da la probabilidad de abandono para cada cliente. Una primera reacción es poner un punto de corte para etiquetar a los clientes como “abandonadores” o “no abandonadores”. Esto no es tan buena idea.

Supongamos que

• el tratamiento de retención cuesta 1200 pesos por cliente,
• estimamos mediante pruebas que nuestro tratamiento reduce la probabilidad de abandono en un 60%
• Tenemos algún tipo de valuación del valor de los clientes.

Usando las probabilidades podemos decidir en estrategias de aplicación del tratamiento. Simulamos una cartera de clientes y sus valuaciones (que suponemos constantes, pero normalmente también son salidas de modelos predictivos). Las probabilidades de abandono suponemos que están dada por un modelo:

# esta tabla nos da la probabilidad de abandono según
# algún modelo que ajustamos
clientes <- tibble(id = 1:20000, valor = 10000) %>% 
    mutate(prob_pred = rbeta(length(valor), 1, 2)) 

calc_perdida <- function(corte, factor_ret, costo){
    perdida_no_trata <- filter(clientes, prob_pred < corte) %>% 
        mutate(costo = ifelse(rbinom(length(prob_pred), 1, prob_pred) == 1, valor, 0)) %>% 
        summarise(total = sum(costo)) %>% 
        pull(total)
    perdida_trata <- filter(clientes, prob_pred >= corte) %>% 
        mutate(costo = ifelse(rbinom(length(prob_pred), 1, prob_pred*factor_ret) == 1, valor, 0)) %>% 
        summarise(total = sum(costo)) %>% 
        pull(total)
    perdida_cf <- filter(clientes, prob_pred >= corte) %>%  
        mutate(costo = ifelse(rbinom(length(prob_pred), 1, prob_pred) == 1, valor, 0)) %>% 
        summarise(total = sum(costo)) %>% 
        pull(total)
    total <- perdida_no_trata +  perdida_trata - (perdida_no_trata + perdida_cf) +
      costo*nrow(filter(clientes, prob_pred > corte)) 
    total
}
perdidas_sim <- map_dfr(rep(seq(0,1, 0.1), 50), 
    function(x){
      perdida_sim <- calc_perdida(x, 0.6, 1000)
      tibble(perdida = perdida_sim, corte = x)
    }) %>% bind_rows 

ggplot(perdidas_sim, aes(x = factor(corte), y = - perdida / 1e6)) + 
  geom_boxplot() + ylab("Ganancia vs Ninguna acción (millones)") +
  xlab("Corte inferior de tratamiento (prob)")

¿Dónde habría que hacer el punto de corte para tratar a los clientes?

5.2 Análisis de error para clasificadores binarios

En muchas ocasiones, los costos y las decisiones todavía no están bien definidas, y requierimos una manera de evaluar los modelos que sea útil para considerar si el desempeño del modelo es apropiado. En estos casos podemos hacer un análisis de error simplificado que ayude a dirigir nuestro trabajo.

Cuando la variable a predecir es binaria (dos clases), podemos etiquetar una clase como positiva y otra como negativa. En el fondo no importa cómo catalogemos cada clase, pero para problemas particulares una asignación puede ser más natural. Por ejemplo, en diagnóstico de enfermedades, positivo=tiene la enfermedad, en análisis de crédito, positivo=cae en impago, en sistemas de recomendacion, positivo = le gusta el producto X, en recuperación de textos, positivo=el documento es relevante a la búsqueda, etc.

Supondremos entonces que hemos construido un clasificador \(\hat{G}_\alpha\) a partir de probabilidades estimadas \(\hat{p}_1(x) > \alpha\). Por ejemplo, podemos construir el clasificador de Bayes clasificando como positivo a todos los casos \(x\) que cumplan \(\hat{p}_1(x) > 0.5\), y negativos al resto.

Hay dos tipos de errores en un clasificador binario (positivo - negativo):

• Falsos positivos (fp): clasificar como positivo a un caso negativo.
• Falsos negativos (fn): clasificar como negativo a un caso positivo.

A los casos clasificados correctamente les llamamos positivos verdaderos (pv) y negativos verdaderos (nv).

La matriz de confusion es entonces

tabla <- tibble(' ' = c('positivo.pred','negativo.pred','total'),
                    'positivo'=c('vp','fn','pos'),
                    'negativo'=c('fp','vn','neg'),
                    'total' = c('pred.pos','pred.neg',''))
knitr::kable(tabla)

	positivo	negativo	total
positivo.pred	vp	fp	pred.pos
negativo.pred	fn	vn	pred.neg
total	pos	neg

Nótese que un clasificador bueno, en general, es uno que tiene la mayor parte de los casos en la diagonal de la matriz de confusión.

Podemos estudiar a nuestro clasificador en términos de las proporciones de casos que caen en cada celda, que dependen del desempeño del clasificador en cuanto a casos positivos y negativos. La nomenclatura puede ser confusa, pues en distintas áreas se usan distintos nombres para estas proporciones:

• Tasa de falsos positivos: proporción de veces que nos equivocamos en la clasificación para casos negativos. \[\frac{fp}{fp+nv}=\frac{fp}{neg}\]

• Tasa de falsos negativos: proporción de veces que nos equivocamos en la clasificación para casos positivos. \[\frac{fn}{pv+fn}=\frac{fn}{pos}\]

• Especificidad: tasa de correctos entre negativos, es el complemento de la tasa de falsos positivos. \[\frac{vn}{fp+vn}=\frac{vn}{neg}\]

• Sensibilidad o Recall o Exhaustividad: tasa de correctos entre positivos, es el complemento de la tasa de falsos negativos. \[\frac{vp}{vp+fn}=\frac{vp}{pos}\]

También existen otras que tienen como base la clasificación en lugar de la clase verdadera, y en algunos casos estas ayudan en la interpretación:

• Precisión o valor predictivo positivo: tasa de correctos entre todos los que clasificamos como positivos \[\frac{vp}{vp+fp}=\frac{vp}{pred.pos}\]

• Valor predictivo negativo: tasa de correctos entre todos los que clasificamos como negativos \[\frac{vn}{fn+vn}=\frac{vn}{pred.neg}\]

Dependiendo de el tema y el objetivo hay medidas más naturales que otras:

• En búsqueda y recuperación de documentos o imagenes, o detección de fraude ( donde positivo = el documento es relevante / la transacción es fraudulenta y negativo = el documento no es relevante / transacción normal), se usa más comunmente precisión y exhaustividad. La exhaustividad (recall o sensibilidad) es muy importante porque mide qué tan completa es la lista de resultados que entregamos. Pero para que el sistema sea utilizable (efectivo para el usuario), la mayor parte de los resultados entregados deben ser relevantes. Entonces es natural usal la precisión, que mide, entre todos los resultados con predicción positiva, que presentamos al usuario, qué porcentaje realmente son relevantes (precisión).

Un clasificador preciso es uno que tal que una fracción alta de sus predicciones positivas son positivos verdaderos. Sin embargo, podría no ser muy sensible: de todos los positivos que hay, solamente clasifica como positivos a una fracción chica. Conversamente, un clasificador podría ser muy sensible: captura una fracción alta de los positivos, pero también clasifica como positivos a muchos casos que son negativos (precisión baja).

• En pruebas para detectar alguna enfermedad, muchas veces se utliza sensibilidad (recall o exhaustividad) y especificidad (qué proporción de negativos excluimos). Sensibilidad nos dice si estamos capturando a la mayoría de los verdaderos positivos. La especificidad nos dice qué tan bien estamos descartando casos que son verdaderos negativos. Si después de la prueba se propone algún tratamiento, estas dos cantidades nos dicen qué proporción de “enfermos” vamos a tratar y qué proporción de “sanos” vamos a tratar.

Cada clasificador tiene un balance distinto precisión y sensibilidad (o sensibilidad y especificidad) . Muchas veces no escogemos clasificadores por la tasa de incorrectos solamente, sino que intentamos buscar un balance adecuado entre el comportamiento de clasificación para positivos y para negativos.

Medidas resumen de desempeño

Medidas resumen populares de los errores de un clasificador particular son:

• Tasa de clasificación incorrecta \[\frac{fn+fv}{neg+pos}\] Y existen otras medidas que intentan resumir los dos tipos de errores de distinta manera, como

• Medida F (media armónica de precisión y recall) \[2\frac{precision \cdot recall}{precision + recall}\] Se usa la la media armónica que penaliza más fuertemente desempeño malo en alguna de nuestras dos medidas (precisión y recall) que el promedio armónico.

Ejemplo (medida F)

Si precision = 0.01 (muy malo) y recall = 1 (excelente), o recall=0.01 y precisión = 1 (excelente), la media usual considera igual de buenos estos dos clasificadores. A su vez, estos dos se califican similar a un clasificador con precision = 0.5 y recall = 0.5. Sin embargo, la media armónica (F) da un score mucho más bajo a los primeros dos clasificadores:

media_armonica <- function(x){
    1/mean(1/x)
}
media_armonica(c(0.01, 1))

## [1] 0.01980198

media_armonica(c(0.5, 0.5))

## [1] 0.5

Nótese que estas medidas no toman en cuenta las probablidades de clase, y son para un clasificador particular.

Ejercicio

Calcular la matriz de confusión (sobre la muestra de prueba) para el clasificador logístico de diabetes en términos de imc y edad (usando punto de corte de 0.5). Calcula adicionalmente con la muestra de prueba sus valores de especificidad y sensibilidad, y precisión y recall.

diabetes_ent <- as_tibble(MASS::Pima.tr) %>% mutate(type = as.character(type))
diabetes_pr <- as_tibble(MASS::Pima.te) %>% mutate(type = as.character(type))
# normalizer
receta_diabetes <- recipe(type ~ ., diabetes_ent) %>%
  step_mutate(type = factor(type, levels = c("Yes", "No"))) %>% 
  step_normalize(all_predictors()) %>% 
  prep()
# ajustar
library(keras)

## 
## Attaching package: 'keras'

## The following object is masked from 'package:yardstick':
## 
##     get_weights

mod_1 <- logistic_reg() %>% 
  set_engine("keras") %>% 
  set_mode("classification") %>%
  set_args(epochs = 250, optimizer = optimizer_sgd(lr = 0.3),
           batch_size = nrow(diabetes_ent), 
           verbose = FALSE) %>% 
  fit(type ~ bmi + age, juice(receta_diabetes))
# otra opcion es
# mod_1 <- logistic_reg() %>% set_engine("glm") %>% set_mode("classification")

Ahora probamos. Primero calculamos directamente:

prueba_baked <- bake(receta_diabetes, diabetes_pr)
preds_prueba <- predict(mod_1, prueba_baked, type ='prob') %>% 
  bind_cols(prueba_baked)
# usar punto de corte 0.5
preds_prueba <- preds_prueba %>% 
  mutate(clase_pred_g = ifelse(preds_prueba$.pred_Yes > 0.5, "pred_Yes", "pred_No")) 
# calcular matriz de confusión
confusion <- preds_prueba %>% 
  group_by(type, clase_pred_g) %>% 
  count() %>% pivot_wider(names_from = type, values_from = n) %>% 
  ungroup() %>% 
  column_to_rownames("clase_pred_g")
# en los renglones están las predicciones
confusion[c("pred_Yes", "pred_No"), ]

##          Yes  No
## pred_Yes  42  29
## pred_No   67 194

Finalmente podemos calcular:

sensibilidad_1 <- confusion["pred_Yes", "Yes"] / sum(confusion[, "Yes"])
precision_1 <- confusion["pred_Yes", "Yes"] / sum(confusion["pred_Yes", ])
sprintf("Precisión: %.2f, Sensibilidad (recall): %.2f", precision_1, sensibilidad_1)

## [1] "Precisión: 0.59, Sensibilidad (recall): 0.39"

O también podemos hacer:

preds_prueba <- preds_prueba %>% 
  mutate(clase_pred = ifelse(.pred_Yes > 0.5, "Yes", "No")) %>% 
  mutate(clase_pred = factor(clase_pred, levels = c("Yes", "No")))
## calcular con yardstick
metricas_1 <- metric_set(accuracy, precision, recall)
preds_prueba %>% 
  metricas_1(type, estimate = clase_pred)

## # A tibble: 3 x 3
##   .metric   .estimator .estimate
##   <chr>     <chr>          <dbl>
## 1 accuracy  binary         0.711
## 2 precision binary         0.592
## 3 recall    binary         0.385

También podemos evaluar con

metricas_2 <- metric_set(accuracy, sens, spec)
preds_prueba %>% 
  metricas_2(type, estimate = clase_pred)

## # A tibble: 3 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.711
## 2 sens     binary         0.385
## 3 spec     binary         0.870

5.3 Análisis de error para probabilidades de clase

Con las métricas interpretables de la sección anterior (sensibilidad, especificidad, etc) es posible construir métricas resumen para evaluar probabilidades de clase. La idea general es:

Si tenemos un modelo que estima probablidades de clase (para un problema binario), podemos construir distintos clasificadores (decisiones) usando distintos puntos de corte. Para cada \(0 < \alpha <1\), tenemos un clasificador \(\hat{G}_{\alpha}\) dado por:

• Si \(p_1(x) > \alpha\), clasificamos el caso como positivo (1). En otro caso clasificamos como negativo.
• Buscamos resúmenes que incorporen desempeño de todos los posibles clasificadores \(\hat{G}_{\alpha}\), usando las medidas de arriba de desempeño para clasificadores binarios.

Distintos valores de de \(\alpha\) dan distintos perfiles de sensibilidad-especificidad para una misma estimación de las probabilidades condicionales de clase: Para minimizar la tasa de incorrectos conviene poner \(\alpha= 0.5\). Sin embargo, es común que este no es el único fin de un clasificador bueno (pensar en ejemplo de fraude).

• Cuando incrementamos \(\alpha\), quiere decir que exigimos estar más seguros de que un caso es positivo para clasificarlo como positivo. Esto quiere decir que la sensibilidad (recall) tiende a ser más chica. Por otro lado la precisión tiende a aumentar, pues el porcentaje de verdaderos positivos entre nuestras predicciones positivas será mayor. También la especificidad tiende a ser más grande.

Ejemplo

Por ejemplo, si en el caso de diabetes incrementamos el punto de corte a 0.7:

preds_prueba <- preds_prueba %>% 
  mutate(clase_pred = ifelse(.pred_Yes > 0.7, "Yes", "No")) %>% 
  mutate(clase_pred = factor(clase_pred, levels = c("Yes", "No")))
## calcular con yardstick
metricas <- metric_set(accuracy, precision, recall)
preds_prueba %>% 
  metricas(type, estimate = clase_pred)

## # A tibble: 3 x 3
##   .metric   .estimator .estimate
##   <chr>     <chr>          <dbl>
## 1 accuracy  binary         0.702
## 2 precision binary         0.656
## 3 recall    binary         0.193

La precisión mejora, pero la sensibilidad (recall) se deteriora. Visto desde el punto de vista de sensiblidad y especificidad:

preds_prueba %>% 
  metricas_2(type, estimate = clase_pred)

## # A tibble: 3 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.702
## 2 sens     binary         0.193
## 3 spec     binary         0.951

• Cuando hacemos más chico \(\alpha\), entonces exigimos estar más seguros de que un caso es negativo para clasificarlo como negativo. Esto aumenta la sensibilidad, pero la precisión baja. Por ejemplo, si en el caso de diabetes ponemos el punto de corte en 0.3:

Ejemplo

preds_prueba <- preds_prueba %>% 
  mutate(clase_pred = ifelse(.pred_Yes > 0.3, "Yes", "No")) %>% 
  mutate(clase_pred = factor(clase_pred, levels = c("Yes", "No")))
## calcular con yardstick
metricas <- metric_set(accuracy, precision, recall)
preds_prueba %>% 
  metricas(type, estimate = clase_pred)

## # A tibble: 3 x 3
##   .metric   .estimator .estimate
##   <chr>     <chr>          <dbl>
## 1 accuracy  binary         0.660
## 2 precision binary         0.488
## 3 recall    binary         0.734

preds_prueba %>% 
  metricas_2(type, estimate = clase_pred)

## # A tibble: 3 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.660
## 2 sens     binary         0.734
## 3 spec     binary         0.623

---

Buscamos ahora hacer análisis de este tipo para todos los posibles puntos de corte. Dos técnicas son curvas ROC y curvas de precisión-exhaustividad, y medidas resumen asociadas.

5.4 Espacio ROC de clasificadores

Podemos visualizar el desempeño de cada uno de estos clasificadores construidos con puntos de corte mapeándolos a las coordenadas de tasa de falsos positivos (1-especificidad) y sensibilidad:

clasif_1 <- tibble(
  corte = c('0.3','0.5','0.7','perfecto','azar'),
  tasa_falsos_pos=c(0.24,0.08,0.02,0,0.7),
  sensibilidad =c(0.66, 0.46,0.19,1,0.7))
ggplot(clasif_1, aes(x=tasa_falsos_pos, y=sensibilidad,
  label=corte)) + geom_point() + 
  geom_abline(intercept=0, slope=1) +
  xlim(c(0,1)) +ylim(c(0,1)) + geom_text(hjust=-0.3, col='red')+
  xlab('1-especificidad (tasa falsos pos)')

1. Nótese que agregamos otros dos clasificadores, uno perfecto, que tiene tasa de falsos positivos igual a 0 y sensibilidad igual a 1.
2. En esta gráfica, un clasificador \(G_2\) que está arriba a la izquierda de \(G_1\) domina a \(G_1\), pues tiene mejor especificidad y mejor sensibilidad. Entre los clasificadores 0.3, 0.5 y 0.7 de la gráfica, no hay ninguno que domine a otro.
3. Todos los clasificadores en la diagonal son equivalentes a un clasificador al azar. ¿Por qué? La razón es que si cada vez que vemos un nuevo caso lo clasificamos como positivo con probabilidad \(p\) fija y arbitraria. Esto implica que cuando veamos un caso positivo, la probabilidad de ’atinarle’ es de p (sensibilidad), y cuando vemos un negativo, la probabilidad de equivocarnos también es de 1-p (tasa de falsos positivos), por lo que la espcificidad es p también. De modo que este clasificador al azar está en la diagonal.
4. ¿Qué podemos decir acerca de clasificadores que caen por debajo de la diagonal? Estos son clasificadores particularmente malos, pues existen clasificadores con mejor especificidad y/o sensibilidad que son clasificadores al azar! Sin embargo, se puede construir un mejor clasificador volteando las predicciones, lo que cambia sensibilidad por tasa de falsos positivos.
5. ¿Cuál de los tres clasificadores es el mejor? En términos de la tasa de incorrectos, el de corte 0.5. Sin embargo, para otros propósitos puede ser razonable escoger alguno de los otros.

5.5 Perfil de un clasificador binario y curvas ROC

En lugar de examinar cada punto de corte por separado, podemos hacer el análisis de todos los posibles puntos de corte mediante la curva ROC (receiver operating characteristic, de ingeniería).

Para un problema de clasificación binaria, dadas estimaciones \(\hat{p}_1(x)\), la curva ROC grafica todos los pares de (1-especificidad, sensibilidad) para cada posible punto de corte \(\hat{p}_1(x) > \alpha\).

Ejemplo

roc_tbl <- roc_curve(preds_prueba, truth = type, estimate = .pred_Yes)
roc_tbl

## # A tibble: 321 x 3
##    .threshold specificity sensitivity
##         <dbl>       <dbl>       <dbl>
##  1  -Inf          0                 1
##  2     0.0525     0                 1
##  3     0.0557     0.00448           1
##  4     0.0602     0.00897           1
##  5     0.0623     0.0135            1
##  6     0.0648     0.0224            1
##  7     0.0655     0.0269            1
##  8     0.0665     0.0314            1
##  9     0.0670     0.0359            1
## 10     0.0700     0.0404            1
## # … with 311 more rows

ggplot(roc_tbl, aes(x = 1 - specificity, y = sensitivity)) +
  geom_path(aes(colour = .threshold), size = 1.2) +
  geom_abline(colour = "gray") + 
  coord_equal() +
  xlab("Tasa de falsos positivos") + ylab("Sensibilidad")

En esta gráfica podemos ver todos los clasificadores posibles basados en las probabilidades de clase. Podemos usar estas curvas como evaluación de nuestros clasificadores, dejando para más tarde la selección del punto de corte, si esto es necesario (por ejemplo, dependiendo de los costos de cada tipo de error).

También podemos definir una medida resumen del desempeño de un clasificador según esta curva:

La medida AUC (area under the curve) para un clasificador es el área bajo la curva generada por los pares sensibilidad-especificidad de la curva ROC.

roc_auc(preds_prueba, type, .pred_Yes)

## # A tibble: 1 x 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 roc_auc binary         0.761

Cuanto más cerca de uno es este valor, mejor discriminan las probabilidades.

También es útil para comparar modelos. Consideremos el modelo de los datos de diabetes que incluyen todas las variables:

mod_2 <- logistic_reg() %>% 
  set_engine("keras") %>% 
  set_mode("classification") %>%
  set_args(epochs = 250, optimizer = optimizer_sgd(lr = 0.3),
           batch_size = nrow(diabetes_ent), 
           verbose = FALSE) %>% 
  fit(type ~ ., juice(receta_diabetes))
preds_prueba_completo <- predict(mod_2, prueba_baked, type ='prob') %>% 
  bind_cols(prueba_baked)
preds_prueba_2 <- bind_rows(preds_prueba %>% mutate(modelo = "IMC y edad"),
                            preds_prueba_completo %>% mutate(modelo = "Todas")) %>% 
  group_by(modelo)

Y graficamos juntas:

roc_2_tbl <- roc_curve(preds_prueba_2, type, .pred_Yes)
ggplot(roc_2_tbl, aes(x = 1 - specificity, y = sensitivity)) +
  geom_path(aes(colour = modelo), size = 1.2) +
  geom_abline(colour = "gray") + 
  coord_equal() +
  xlab("Tasa de falsos positivos") + ylab("Sensibilidad")

Comparación auc:

roc_auc(preds_prueba_2, type, .pred_Yes) 

## # A tibble: 2 x 4
##   modelo     .metric .estimator .estimate
##   <chr>      <chr>   <chr>          <dbl>
## 1 IMC y edad roc_auc binary         0.761
## 2 Todas      roc_auc binary         0.866

En este ejemplo, vemos que casi no importa qué perfil de especificidad y sensibilidad busquemos: el clasificador que usa todas las variables domina casi siempre al clasificador que sólo utiliza IMC y edad.

5.6 Curvas de precisión-sensibilidad

Para mostrar los posibles perfiles de clasificación de nuestro modelo, podemos mostrar las posibles combinaciones de precisión y recall bajo todos los posibles cortes:

pr_auc_tbl <- pr_curve(preds_prueba_2, type, .pred_Yes)
ggplot(pr_auc_tbl %>% filter(modelo != "Todas"), 
       aes(x = recall , y = precision)) +
  geom_path(aes(colour = .threshold), size = 1.2) +
  geom_abline(colour = "gray") + 
  coord_equal() +
  xlab("Exhaustividad (Sensibilidad)") + ylab("Precisión")

Nota que hay algunas oscilaciones en la curva. Cuando el punto de corte sube, la sensibilidad siempre baja o se queda igual (hay la misma cantidad o menos de verdaderos positivos). Sin embargo, la precisión se calcula sobre la base del número predicciones positivas, y esta cantidad siempre disminuye cuando el punto de corte aumenta. Especialmente cuando el número de predicciones positivas es chico, esto puede producir oscilaciones como las de la figura.

Ahora probamos usando todas las variables y comparamos

ggplot(pr_auc_tbl , 
       aes(x = recall , y = precision)) +
  geom_path(aes(colour = modelo), size = 1.2) +
  xlab("Exhaustividad (Sensibilidad)") + ylab("Precisión")

Y observamos que el modelo con todas las variables tiende a domina al modelo de IMC y edad.

Observaciones: - Precisión y Recall (sensibilidad) son medidas muy naturales para muchos problemas, en particular aquellos donde la clase de positivos es relativamente chica. - Sin embargo, en puntos de corte altos las estimaciones de precisión son ruidosas cuando nuestro conjunto de prueba es relativamente chico y no existen muchos casos con probabilidades altas (el denominador de precisión es el número de clasificados como positivos) - Una alternativa es usar curvas ROC

5.7 Curvas acumuladas de ganancia y lift

Otro enfoque común y útil es graficar en el eje horizontal el porcentaje de casos que cada clasificador clasifica como positivos para cada punto de corte- en términos de la decisión que nos interesa, esto usualmente indica qué porcentaje de los casos esperamos intervenir o probar (por ejemplo, el porcentaje que recibe una campaña publicitaria, o que califica para algún programa especial). En el eje vertical ponemos la sensibilidad, que es qué porcentaje de todos los positivos esperamos encontrar:

gain_tbl <- gain_curve(preds_prueba_2, type, .pred_Yes)
ggplot(gain_tbl, 
       aes(x = .percent_tested , y = .percent_found)) +
  geom_path(aes(colour = modelo), size = 1.2) +
  geom_abline(colour = "gray") + 
  coord_equal() +
  xlab("% Clasificado positivo") + ylab("Sensibilidad")

Y un análogo para las curvas de precisión-recall es la curva del lift. En el eje vertical graficamos la sensibilidad entre el % de clasificados como positivos (que es lo mismo que la precisión entre el % de casos positivos):

lift_tbl <- lift_curve(preds_prueba_2, type, .pred_Yes)
ggplot(lift_tbl, 
       aes(x = .percent_tested , y = .lift)) +
  geom_path(aes(colour = modelo), size = 1.2) +
  xlab("% Clasificado positivo") + ylab("Lift")

## Warning: Removed 2 row(s) containing missing values (geom_path).

Cuando el lift es más grande que uno, quiere decir que el modelo es superior en encontra positivos que una estrategia de tomar casos al azar. Por las mismas razones que las curvas de precisión-exhaustividad, la gráfica puede ser ruidosa para puntos de corte altos.

Regresión logística para más de 2 clases

Consideramos ahora un problema con más de dos clases, de manera que \(G ∈ {1,2,...,K}\) (\(K\) clases), y tenemos \(X = (X1 ...,Xp)\) entradas. ¿Cómo generalizar el modelo de regresión logística para este problema? Una estrategia es la de uno contra todos:

En clasificación uno contra todos, hacemos

1. Para cada clase \(g\in\{1,\ldots,K\}\) entrenamos un modelo de regresión logística (binaria) \(\hat{p}^{(g)}(x)\), tomando como positivos a los casos de 1 clase \(g\), y como negativos a todo el resto. Esto lo hacemos como en las secciones anteriores, y de manera independiente para cada clase.

2. Para clasificar un nuevo caso \(x\), calculamos \[\hat{p}^{(1)}, \hat{p}^{(2)},\ldots, \hat{p}^{(K)}\]

y clasificamos a la clase de máxima probabilidad \[\hat{G}(x) = \arg\max_g \hat{p}^{(g)}(x)\] Nótese que no hay ninguna garantía de que las probabilidades de clase sumen 1, pues se trata de estimaciones independientes de cada clase. En este sentido, produce estimaciones que en realidad no satisfacen las propiedades del modelo de probabilidad establecido — aunque pueden normalizarse. Sin embargo, esta estrategia es simple y en muchos casos funciona bien.

5.8 Regresión logística multinomial

Si queremos obtener estimaciones de las probabilidades de clase que sumen uno, entonces tenemos que contruir las estimaciones de cada clase de clase de manera conjunta. Como vimos antes, tenemos que estimar, para cada \(x\) y \(g\in\{1,\ldots, K\}\), las probabilidades condicionales de clase: \[p_g(x) = P(G = g|X = x).\] Una manera simple de construir estas clases es considerando un juego de coeficientes para cada clase, por ejemplo, para la clase \(j\) consideramos el predictor lineal dado por:

\[L_j = \beta_0^j + \beta_1^jx_1 + \ldots + \beta_p^j x_p.\] No podemos sumar y normalizar estos valores pues pueden tomar valores negativos, por ejemplo. Esto podemos evitarlo tomando la exponencial:

\[Z_j = \exp(\beta_0^j + \beta_1^jx_1 + \ldots + \beta_p^j x_p).\] Estos valores \(Z_j\) son positivos, así que podemos definir

\[Z = \sum_{j=1}^L Z_j, \] y calcular

\[p_j(x) = \frac{Z_j}{Z} = \frac{\exp(\beta_0^j + \beta_1^jx_1 + \ldots + \beta_p^j x_p)}{Z}\] A estas probabilidades (que suman 1 sobre todas las clases) se les llama el softmax de las cantidades \(L_j\), es decir, el softmax de \(L_1, \ldots, L_K\) es

\[\frac{\exp(L_j)}{\sum_{j=1}^K \exp(L_j)}\] Esta es la idea básica de regresión logística multinomial.

Nótese sin embargo que estos coeficientes sobreparametrizan las probabilidades, pues como estas probabilidades suman 1, una vez que conocemos las primeras \(K-1\), la otra debería estar determinada. Esto implica, entre otras cosas, que hay múltiples formas de expresar con el mismo modelo con distintos parámetros (¿ejemplo?).

• Este es un problema que vamos a resolver más adelante usando regularización. Esta sobreparametrización se usa frecuentemente en modelos de redes neuronales, por ejemplo.
• Sin embargo, también es posible escoger una clase de referencia (por ejemplo la \(K\)), y poner \(\beta^K_{j}=0\) para toda \(j\), de forma que \(Z_K=1\), y entonces

\[p_K(x) = \frac{1}{1 + \sum_{j=1}^{K-1} Z_j},\] y las probabilidades suman uno otra vez. Esta parametrización es más frecuente en regresión logística multinomial tradicional.

• En regresión logística usamos esta última parametrización: para dos clases, tenemos que \[p_1(x) = \frac{\exp(\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_o)}{1 + \exp(\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p)} \] y como \(K=2\), \[p_2(x) = \frac{1}{1 + \exp(\beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p)}\]

Una vez definido el modelo multinomial, podemos usar algún método númerico para minimizar la devianza (por ejemplo, descenso en gradiente).

Interpretación de coeficientes

Para comparar la clase \(j\) con la clase \(k\) notamos que

\[\log\left (\frac{p_j(x)}{p_k(x)}\right ) = (\beta_{0,{j}}- \beta_{0,{k}}) + (\beta_{1,{j}}-\beta_{1,{k}} )x_1 + \ldots + (\beta_{p,{j}} -\beta_{p,{k}}) x_p\]

Así que sólo hace falta restar los coeficientes. Estos coeficientes nos dicen cómo cambia \(p_j\) en relación a \(p_k\) cuando una entrada \(x_i\) cambia.

5.9 Ejemplo: Clasificación de dígitos con regresión multinomial

digitos_entrena <- read_csv('../datos/zip-train.csv')
digitos_prueba <- read_csv('../datos/zip-test.csv') 
names(digitos_entrena)[1] <- 'digito'
names(digitos_entrena)[2:257] <- paste0('pixel_', 1:256)
names(digitos_prueba)[1] <- 'digito'
names(digitos_prueba)[2:257] <- paste0('pixel_', 1:256)
digitos_entrena <- digitos_entrena %>% 
  mutate(digito = factor(digito, levels = seq(0, 9, 1)))
digitos_prueba <- digitos_prueba %>% 
  mutate(digito = factor(digito, levels = seq(0, 9, 1)))

En este ejemplo, usamos la función multinom de nnet, que usa BFGS para hacer la optimización (puedes también usar keras o glmnet que veremos más adelante):

reg_multinomial <- multinom_reg() %>% 
  set_engine("nnet") %>% 
  set_args(MaxNWt = 100000, maxit = 20)
mod_mult <- reg_multinomial %>% fit(digito ~ ., data = digitos_entrena)

Checamos para diagnóstico la matriz de confusión y devianza de entrenamiento

preds_entrena <- predict(mod_mult, digitos_entrena) %>% 
  bind_cols(digitos_entrena %>% select(digito))
preds_entrena_proba <- predict(mod_mult, digitos_entrena, type = "prob") %>% 
  bind_cols(digitos_entrena %>% select(digito))
preds_entrena_proba %>% mn_log_loss(digito, .pred_0:.pred_9)

## # A tibble: 1 x 3
##   .metric     .estimator .estimate
##   <chr>       <chr>          <dbl>
## 1 mn_log_loss multiclass     0.205

table(preds_entrena$.pred_class, preds_entrena$digito)

##    
##        0    1    2    3    4    5    6    7    8    9
##   0 1153    0    5    2    3    9    1    1    7    0
##   1    0  998    0    0    2    0    1    1    2    3
##   2    2    0  693    1    7    2    8    3   10    2
##   3    9    0   15  632    2   21    0    2   24    2
##   4    3    2    9    2  621    4    4    9   10   44
##   5   24    4    5   19    9  511   43    1   34    6
##   6    2    0    0    0    1    3  607    0    0    0
##   7    0    0    1    0    0    1    0  613    1    8
##   8    1    1    3    2    2    4    0    1  451    2
##   9    0    0    0    0    5    1    0   14    3  577

Ahora validamos con la muestra de prueba y calculamos error de clasificación:

preds <- predict(mod_mult, digitos_prueba) %>% 
  bind_cols(digitos_prueba %>% select(digito))
preds_proba <- predict(mod_mult, digitos_prueba, type = "prob") %>% 
  bind_cols(digitos_prueba %>% select(digito))
preds_proba %>% mn_log_loss(digito, .pred_0:.pred_9)

## # A tibble: 1 x 3
##   .metric     .estimator .estimate
##   <chr>       <chr>          <dbl>
## 1 mn_log_loss multiclass     0.649

confusion_prueba <- table(preds$.pred_class, preds$digito)
confusion_prueba

##    
##       0   1   2   3   4   5   6   7   8   9
##   0 335   0   3   0   3   6   4   0   3   0
##   1   0 252   0   0   1   0   0   0   1   3
##   2   1   1 171   4   8   0   5   2   3   2
##   3   3   3   8 145   1  18   0   3  12   0
##   4   3   6   7   1 176   1   2   6   8  16
##   5  11   1   5  13   2 130  13   2  14   1
##   6   4   1   1   0   2   0 143   0   1   0
##   7   0   0   1   1   2   1   0 130   1   5
##   8   1   0   2   1   2   1   3   0 118   0
##   9   1   0   0   1   3   3   0   4   5 150

sum(diag(confusion_prueba))/sum(confusion_prueba)

## [1] 0.8719482

round(prop.table(confusion_prueba, 2), 2)

##    
##        0    1    2    3    4    5    6    7    8    9
##   0 0.93 0.00 0.02 0.00 0.01 0.04 0.02 0.00 0.02 0.00
##   1 0.00 0.95 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.02
##   2 0.00 0.00 0.86 0.02 0.04 0.00 0.03 0.01 0.02 0.01
##   3 0.01 0.01 0.04 0.87 0.00 0.11 0.00 0.02 0.07 0.00
##   4 0.01 0.02 0.04 0.01 0.88 0.01 0.01 0.04 0.05 0.09
##   5 0.03 0.00 0.03 0.08 0.01 0.81 0.08 0.01 0.08 0.01
##   6 0.01 0.00 0.01 0.00 0.01 0.00 0.84 0.00 0.01 0.00
##   7 0.00 0.00 0.01 0.01 0.01 0.01 0.00 0.88 0.01 0.03
##   8 0.00 0.00 0.01 0.01 0.01 0.01 0.02 0.00 0.71 0.00
##   9 0.00 0.00 0.00 0.01 0.01 0.02 0.00 0.03 0.03 0.85

El resultado no es muy bueno. Veremos más adelante mejores métodos para este problema. ¿Podemos interpretar el modelo?

Una idea es tomar los coeficientes y graficarlos según la estructura de las imágenes:

coefs <- coef(mod_mult$fit)
coefs_reng <- coefs[1, , drop =FALSE]
coefs <- rbind(coefs_reng, coefs)
coefs[1 , ] <- 0
dim(coefs)

## [1]  10 257

beta_df <- coefs[,-1] %>% as.data.frame %>% 
  mutate(digito = 0:(nrow(coefs)-1)) %>%
  gather(pixel, valor, contains('pixel')) %>%
  separate(pixel, into = c('str','pixel_no'), sep='_') %>%
  mutate(x = (as.integer(pixel_no)-1) %% 16, y = -((as.integer(pixel_no)-1) %/% 16))
head(beta_df)

##   digito   str pixel_no        valor x y
## 1      0 pixel        1  0.000000000 0 0
## 2      1 pixel        1  0.621681333 0 0
## 3      2 pixel        1 -0.005914605 0 0
## 4      3 pixel        1  0.044257959 0 0
## 5      4 pixel        1  0.190966643 0 0
## 6      5 pixel        1 -0.010655932 0 0

Podemos cruzar la tabla con sí misma para hacer comparaciones de cómo discrimina el modelo entre cada par de dígitos:

tab_coef <- beta_df %>% select(digito, x, y, valor)
tab_coef_1 <- tab_coef
names(tab_coef_1) <- c('digito_1','x','y','valor_1')
tab_cruzada <- full_join(tab_coef_1, tab_coef) %>% mutate(dif = valor_1 - valor)

## Joining, by = c("x", "y")

tab_cruzada <- tab_cruzada %>% group_by(digito, digito_1) %>% 
  mutate(dif_s = (dif - mean(dif))/sd(dif)) %>%
  mutate(dif_p = pmin(pmax(dif_s, -2), 2))

ggplot(tab_cruzada, aes(x=x, y=y)) + geom_tile(aes(fill = dif_p)) + 
  facet_grid(digito_1~digito) + scale_fill_distiller(palette = "Spectral")

Discusión

Nótese que no corrimos el modelo hasta convergencia. Vamos a hacerlo ahora:

mod_mult <- reg_multinomial %>% 
  set_args(maxit = 500) %>% 
  fit(digito ~ ., data = digitos_entrena)

preds_entrena <- predict(mod_mult, digitos_entrena) %>% 
  bind_cols(digitos_entrena %>% select(digito))
preds_entrena_proba <- predict(mod_mult, digitos_entrena, type = "prob") %>% 
  bind_cols(digitos_entrena %>% select(digito))
preds_entrena_proba %>% mn_log_loss(digito, .pred_0:.pred_9)

## # A tibble: 1 x 3
##   .metric     .estimator    .estimate
##   <chr>       <chr>             <dbl>
## 1 mn_log_loss multiclass 0.0000000116

table(preds_entrena$.pred_class, preds_entrena$digito)

##    
##        0    1    2    3    4    5    6    7    8    9
##   0 1194    0    0    0    0    0    0    0    0    0
##   1    0 1005    0    0    0    0    0    0    0    0
##   2    0    0  731    0    0    0    0    0    0    0
##   3    0    0    0  658    0    0    0    0    0    0
##   4    0    0    0    0  652    0    0    0    0    0
##   5    0    0    0    0    0  556    0    0    0    0
##   6    0    0    0    0    0    0  664    0    0    0
##   7    0    0    0    0    0    0    0  645    0    0
##   8    0    0    0    0    0    0    0    0  542    0
##   9    0    0    0    0    0    0    0    0    0  644

Ahora validamos con la muestra de prueba y calculamos error de clasificación:

preds <- predict(mod_mult, digitos_prueba) %>% 
  bind_cols(digitos_prueba %>% select(digito))
preds_proba_sobre <- predict(mod_mult, digitos_prueba, type = "prob") %>% 
  bind_cols(digitos_prueba %>% select(digito))
preds_proba_sobre %>% mn_log_loss(digito, .pred_0:.pred_9)

## # A tibble: 1 x 3
##   .metric     .estimator .estimate
##   <chr>       <chr>          <dbl>
## 1 mn_log_loss multiclass      5.21

confusion_prueba <- table(preds$.pred_class, preds$digito)
confusion_prueba

##    
##       0   1   2   3   4   5   6   7   8   9
##   0 332   0   6   2   4   2   1   2   7   2
##   1   0 242   1   3   3   4   1   2   0   2
##   2   2   2 148   5   5   0   4   3   3   0
##   3   4   1   9 128   4  10   0   3   2   4
##   4   3   5   8   0 149   8   6   7   5   2
##   5   0   1   3  11   5 116   8   0  10   1
##   6   5   7   4   3  10   4 144   0   4   1
##   7   2   1   3   1   4   1   1 125   2   4
##   8   6   3  14   7   6  10   4   0 132   3
##   9   5   2   2   6  10   5   1   5   1 158

sum(diag(confusion_prueba))/sum(confusion_prueba)

## [1] 0.8340807

round(prop.table(confusion_prueba, 2), 2)

##    
##        0    1    2    3    4    5    6    7    8    9
##   0 0.92 0.00 0.03 0.01 0.02 0.01 0.01 0.01 0.04 0.01
##   1 0.00 0.92 0.01 0.02 0.01 0.03 0.01 0.01 0.00 0.01
##   2 0.01 0.01 0.75 0.03 0.03 0.00 0.02 0.02 0.02 0.00
##   3 0.01 0.00 0.05 0.77 0.02 0.06 0.00 0.02 0.01 0.02
##   4 0.01 0.02 0.04 0.00 0.74 0.05 0.04 0.05 0.03 0.01
##   5 0.00 0.00 0.02 0.07 0.03 0.72 0.05 0.00 0.06 0.01
##   6 0.01 0.03 0.02 0.02 0.05 0.03 0.85 0.00 0.02 0.01
##   7 0.01 0.00 0.02 0.01 0.02 0.01 0.01 0.85 0.01 0.02
##   8 0.02 0.01 0.07 0.04 0.03 0.06 0.02 0.00 0.80 0.02
##   9 0.01 0.01 0.01 0.04 0.05 0.03 0.01 0.03 0.01 0.89

Y nota que la devianza y el error de clasificación son más grande que cuando nos detuvimos con menos iteraciones. ¿Por qué pasa esto? Discute en clase:

• Grafica los coeficientes para este segundo modelo
• ¿En cuál de los dos modelos es más fácil interpretar los coeficientes? ¿En cuál es menor el error?
• ¿Cuál crees que es el problema de este segundo modelo comparado con el primero? ¿Por qué crees que sucede? ¿Cómo podríamos corregir este problema?

5.10 Curvas ROC multiclase

En un problema multiclase, podemos construir curvas individuales para cada una de clases de diferentes maneras. Una primera técnica es considerar la predicción de cada clase contra el resto_

curvas_digitos <- roc_curve(preds_proba_sobre, digito, .pred_0:.pred_9)
curvas_digitos

## # A tibble: 16,679 x 4
##    .level   .threshold specificity sensitivity
##    <chr>         <dbl>       <dbl>       <dbl>
##  1 0      -Inf              0                1
##  2 0         0.             0                1
##  3 0       Inf.Nae-324      0.0419           1
##  4 0         2.80e-322      0.0425           1
##  5 0         6.90e-322      0.0431           1
##  6 0         2.02e-321      0.0437           1
##  7 0         1.82e-319      0.0443           1
##  8 0         5.64e-317      0.0449           1
##  9 0         1.31e-315      0.0455           1
## 10 0         4.20e-314      0.0461           1
## # … with 16,669 more rows

Que podemos graficar como:

autoplot(curvas_digitos)

Donde vemos que las clases más fáciles son 0, 1, 6 y 9. El resto de las clases es más difícil de discriminar. Podemos usar el promedio de estas AUC para dar un resumen de desempeño (ya sea ponderado o no por la frecuencia de cada clase):

# Primer modelo
roc_auc(preds_proba, digito, .pred_0:.pred_9, estimator = "macro") %>% 
  mutate(across(where(is.numeric), round, 3))

## # A tibble: 1 x 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 roc_auc macro          0.986

# Segundo modelo
roc_auc(preds_proba_sobre, digito, .pred_0:.pred_9, estimator = "macro")  %>% 
  mutate(across(where(is.numeric), round, 3))

## # A tibble: 1 x 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 roc_auc macro           0.96