Algunas reflexiones sobre el conjunto de datos Iris y el análisis discriminante

Hace ya la friolera de casi 90 años, en 1936 R. A. Fisher publica el artículo “The Use of Multiple Measurements in Taxonomic Problems”.

Aunque el conjunto de datos ha sido extremadamente explotado a lo largo de la historia de la estadística, una de las cosas que resalto es que cuando se lee el paper original, hay un diagrama que me encanta y que es el siguiente:

Es un ejemplo que pongo en alguna de las clases que imparto en el Master Executive en Advanced Analytic & Data Sciencie de la UIC en https://www.uic.es/en/estudis-uic/business-administration/masters-degree-data-science-advanced-analytics

Es un dibujo ¡Hecho a mano! porque claro, en 1936 los ordenadores no existían y ni se le esperaban, estábamos en una situación pre-bélica y el mundo estaba más pendiente de lo que pasaba con Hitler que en pensar en la IA por ejemplo.

Lo que me sorprende de todo esto, aparte de la técnica que comento después, es el proceso de recolección de la información. Cualquiera que vea el dataset, que obtiene rápidamente a partir de R o Python, observa que se generan unos 125 registros de medidas tomadas sobre 125 flores de Iris que Fisher seleccionó. Me gusta hacer especial hincapié en ¿Cómo lo hizo Fisher? Yo me imagino a un señor bien vestido de la época que se iba al campo con un cuadernito, una lupa y algo para medir e iba buscando distintas variedades de flores hasta completar una tabla tal como la siguiente que también se puede observar en el paper que he mencionado:

Esta tabulación, obviamente viene ya por defecto tanto en R como en Python y la usamos sin darle mayor importancia, como si tuviera que estar ahí porque sí y siendo este por ejemplo el aspecto de algunos de los registros que están en R:

Pues bien, Fisher no sólo tuvo la "santa paciencia" de recoger y anotar los datos "en su libretilla", sino que además en ese artículo que trataba de taxonomía, desarrolla el análisis discriminante para cuya introducción en las clase que lo trato recurro a una referencia mucho más moderna como es An introduction to statistical learning del 2013, donde con notación más clara se desarrolla y explica de un modo muy sencillo como Fischer llega a establecer la siguiente regla de discriminación para el caso de 2 poblaciones:

Cuando explico esta fórmula, me gusta dar contexto a los datos y para ello me imagino que soy Fisher, pero con el inconveniente de que no estoy en el campo, ni tengo el tiempo para ir recogiendo y midiendo flores por ahí, son muy de ciudad, ... Por lo que tengo que simular datos, en este sentido desarrollo un programa en R que consta de 2 partes:

• Una parte desconocida para el investigador y que es el proceso de generación de la información, es decir la simulación que se puede hacer del siguiente modo:

  
  

  #S simulan dos conjuntos de 1000 datos acordes a dos distribuciones

  #normales multivariantes

  
  

  n <- 2000

  p <- 2

  
  

  #Parámetros estadísticos de los datos

  sigmapos <- matrix(c(0.5,0.4,0.4,0.5), 2)

  sigmaneg <- matrix(c(0.5,0.2,0.2,0.5), 2)

  meanpos <- c(1,1)

  meanneg <- c(2,1)

  
  

  #Número de datos positivos y negativos

  
  

  npos <- round(n / 2)

  nneg <- n - npos

  
  

  #Generación de datos

  library(MASS)

  set.seed(1)

  xpos <- mvrnorm(npos, meanpos, sigmapos)

  xneg <- mvrnorm(nneg, meanneg, sigmaneg)

  x <- rbind(xpos, xneg)

  y <- matrix(c(rep(1, npos), rep(-1,nneg)))

  
  

  plot(x, col=ifelse(y>0, 1, 2))

Que genera el siguiente gráfico:

• Y ahora que se tienen "datos para investigar" hay que comprobar si con la fórmula dada por Fisher se puede o no separar las poblaciones que visualmente se ve que existiría una línea (plano discriminante), pero que debe ser encontrada. Por tanto lo que resalto es que el investigador lo que realmente conoce son estos datos finales (aka simulados) y debe calcular los parámetros de la fórmula por tanto mediante código vamos a calcular las medias observadas y las covarianzas observadas para obtener las matrices correspondientes que son bien distintas de las reales y que en el anterior código se han denominado meanpos, meanneg, sigmapos, sigmaneg. Eso sí previamente separamos en training - test porque en una investigación, no se tendría el conocimiento de si lo hecho en una muestra tendría una mínima base de extrapolación:

  
  

  #1.- Construcción Training - Validation

  
  

  ntrain <- round(n * 0.7)

  tindex <- sample(n, ntrain)

  xtrain <- x[tindex,]

  xtest <- x[-tindex,]

  ytrain <- y[tindex]

  ytest <- y[-tindex]

  
  

  training <- data.frame(cbind(xtrain, ytrain))

  names(training)<-c("x1", "x2", "y")

  validation <- data.frame(cbind(xtest, ytest))

  names(validation)<-c("x1", "x2", "y")

  
  

  xytrain <- cbind(xtrain, ytrain)

  xtrainpos <- xytrain[which(ytrain ==1),][,c(1,2)]

  ytrainpos <- xytrain[which(ytrain ==1),][,c(3)]

  xtrainneg<- xytrain[which(ytrain ==-1),][,c(1,2)]

  ytrainneg<- xytrain[which(ytrain ==-1),][,c(2)]

  
  

  #¿Qué ocurriría si se reestiman los vectores de medias de las distribuciones con

  #estas simulaciones con los datos de varianza-covarianzas?

  
  

  media1 = c(mean(xtrainpos[,1]),mean(xtrainpos[,2]))

  media2 = c(mean(xtrainneg[,1]),mean(xtrainneg[,2]))

  
  

  desviacion1 <- cov(xtrainpos)

  desviacion2 <- cov(xtrainneg)

  
  

  prior1 = length(ytrainpos) / (length(ytrainpos) + length(ytrainneg))

  prior2 = length(ytrainneg) / (length(ytrainpos) + length(ytrainneg))

  
  

Lo que realmente conozco es por tanto los vectores media1, media2 y la matrices de varianzas-covarianzas observada de cada población desviacioin1 y desviacion2 y que supongo representantes de la población a estudiar ya que recuerdo, en una situación real, no conocería en general los vectores de medias y las matrices de varianzas-covarianzas reales (aka simuladas).

Aunque en R el uso de una librería como le e1071 te permitiría el cálculo de un linear discriminant análisis (que era la propuesta simplificada de Fisher ya que en la época la fórmula que él dio resultaba aún complicada de manejar, aunque R dispone también de la qda para la fórmula al completo) según:

library(e1071)

Modelo_LDA <- lda(y~x1 + x2, data = training, prior =c(0.5, 0.5))

Prediction_LDA <- predict(Modelo_LDA, validation)$class

#Confusion Matrix

table(`Actual Class` = ytest, `Predicted Class` = Prediction_LDA)

error.rate.Prediction_LDA <- sum(ytest != Prediction_LDA)/length(ytest)

print(paste0("Accuary (Precision): ", 1 - error.rate.Prediction_LDA))

Yo prefiero plantear previamente a lo anterior la función según su estado original donde se observa, cómo se incorporan todos los elementos:

Fischer <- function(x, mu1, mu2, sigma1, sigma2, pr1, pr2){

LD <- sign(-0.5* (x - mu1) %*% solve(sigma1) %*% (x - mu1) +

0.5* (x - mu2) %*% solve(sigma2) %*% (x - mu2) +

0.5* log(det(sigma2)/det(sigma1), base = exp(1)) -

log(pr2/pr1, base = exp(1)))

return(LD)

}

FischerPred <- vector("numeric", length(ytest))

#Aplicación en test

for (i in 1:length(ytest))

{

FischerPred[i] <- Fischer(xtest[i,], media1, media2, desviacion1, desviacion2, prior1, prior2)

}

#Confusion Matrix

table(`Actual Class` = ytest, `Predicted Class` = FischerPred)

error.rate.FischerPred <- sum(ytest != FischerPred)/length(ytest)

print(paste0("Accuary (Precision): ", 1 - error.rate.FischerPred))

Aunque si se observan las matrices de confusión, existen algunas diferencias, estas son mínimas y están provocadas por la elección de las prior que "si fuésemos por el camino fácil", al no tener mucho conocimiento, los valores 0.5 y 0.5 se muestran razonables.

Por tanto y para finalizar, con esta clase lo que pretendo enseñar es:

1. El alumno tome conciencia de que los datos vienen de algún sitio y antes de tabularlos, no siempre son fáciles de conseguir

   
   

2. Que se entienda que por un lado está la población real, que será desconocida y por otro está lo que observamos, una muestra que fluctúa y que se espera, represente a la población real, claramente si Fisher hubiera tomado 150 flores el resultado hubiera sido análogo, pero los número distintos, no fijos

   
   

3. También que las fórmulas y/o modelos que están en R tienen un sustrato matemático que ya nos lo sirve el programa de un modo sencillo, pero si no se entiende lo que hay por detrás, estamos usando cajas totalmente negras sin entender el por qué

Espero que estas notas sean de utilidad a alguien y si habéis llegado hasta aquí, al menos os haya motivado a aprender más sobre análisis discriminante que aunque muchos puedan pensar que es bastante viejuno, ¡Oh sorpresa! Son la base de los modelos de Altman de finales de los 60 básicos para entender el riesgo business to business y por supuesto es otra técnica más para complementar aquellas que tratan de clasificación o de separación de poblaciones como la regresión logística, decision trees, ...

LO QUE NO TE CUENTA TU CONSULTOR: Tú eres responsable de tu modelo

Desde que me empecé a dedicar al riesgo de crédito allá por el 2009 y tras estar unos 8-9 años imbuido en el mundo de las series temporales, he estado en diversas empresas tanto en cliente final como en proyectos bajo distintas consultoras (y sí, yo también he sido partícipe de la venta de humo de distintas calidades ...) y lo que me encuentro cuando llego a un proyecto de modelización en una empresa usuaria (o cliente) es generalmente el siguiente panorama:

• Un proveedor externo ha desarrollado, e incluso ejecuta un modelo que se usa en producción (cuando el externo ejecuta y cobra por ejecución del modelo, es lo el éxtasis de la consultoría porque entonces el cliente es lo que se denomina, una "vaca lechera")

  
  

• Apoyados en ocasiones por mi queridísima GDPR (nótese el sarcasmo ...) y demás eufemismos como (el know-how), ellos mismos, constructores del modelo, le hacen un seguimiento mediocre que, ¡Oh, sorpresa! ¡Sale estupendo! en algún parámetro que medianamente medio entiende la alta dirección como el Gini o la RoC (podéis preguntarle a los directivos aprueban las implantaciones que os lo describan con detalle y os echáis unas risas)

  
  

• En general, la empresa cliente se desentiende y no le dedica apenas tiempo, ni tan siquiera a corroborar lo más básico hasta que hay un momento de reducción de costes o se ha liado algo importante y entonces a lo más que se llega, de modo reactivo es a decir que sale muy caro el modelo, pero tras una inversión más o menos cuantiosa y tras unos 4 o 5 años ... ¡No hay narices a quitarlo!

Pues bien lo primero que hago y la vedad, que dada vez me aburre más porque pasa el tiempo y no se aprende, es montar el puñetero sistema de monitoring que debería estar montado desde el momento inicial y esto no es responsabilidad del consultor, esto es responsabilidad de la empresa cliente que ingenuamente se cree lo que le dice el que el vendor model dice: "te vamos a ahorrar trabajo", "te vamos a ahorrar costes", ... ya sabéis esas palabras mágicas que todo directivo - CEO le encantan.

Por tanto voy a tratar un poco qué hacer en el caso de Riesgos ya que es el que más me toca por ahora y desde luego, aunque en la gran banca está generalmente resuelto, fuera de ahí, la fiesta es divertida pero en todo caso, es el área donde más atención se pone, porque si comentamos sobre las áreas Comerciales donde están los modelos de Onboarding, Contratación, Vinculación, Fuga, ... eso sí que que desde luego es un auténtico desastre y que ya veré si tengo energía de hacer un coleccionable sobre sus miserias porque da para varias sesiones.

Lo primero que uno tiene que hacer cuando llega, si se va a dedicar al tema de gestionar modelos es hacer la pregunta ¿Cuántos modelos tenemos y dónde está la documentación? Depende de a quién haga la pregunta, puede desde no encontrar respuesta hasta encontrar distintas respuesta contradictorias y con un poco de suerte, en algún caso le darán lo que alguien entiende por documentación que será o un word o power point mal hecho de unas pocas hojas, por lo que ya desde ese momento, nota que la cosa va ser difícil.

Por tanto hay un trabajo duro en encontrar por cada modelo en responder lo siguiente:

• ¿Dónde están las puntuaciones (scoring, probabilidades, ...) de los modelos?

• ¿Cómo asociarla a las poblaciones?

• Y lo más complejo: como re-construir retrospectivamente la target

Con lo anterior se empezaría a solucionar sólo una parte del seguimiento de modelos y de lo que se trataría sería de conseguir unas gráficas de este estilo (con datos totalmente inventados por mí) que en cierto modo te genera tu vendor model, para venderte más en un futuro pero que ¡DEBERÍAS HACERLAS TÚ!:

Pues bien, lo anterior hacerlo por cada uno de los modelos que tengáis ya os digo que es un reto que hay que conseguir visualizar y mostrar y explicar al resto de responsables - directivos para tomar conciencia, pero además, si queréis analizar el nivel variable a través de Information Values, PSI, etc; la dificultad crece exponencialmente.

Tres son los elementos que se requiere para montar un sistema:

1. Lo primero mucha paciencia y ganas de hablar con mucha gente hasta que se encuentre lo poco que hay de documentación o en ocasiones de restos de código, vosotros tendréis que completar lo que falte

   
   

2. Después montar queries para sacar, en el caso más sencillo (sin variables), los gráficos anteriores, 3 son las columnas necesarias un identificador de elemento de la población, una puntuación, y el valor binario de la target a ser ordenada

   
   

3. Salvo que se tenga aplicativo (que generalmente suelen ser caros o malos), con R o con Python va a ser más que suficiente. Yo recomiendo 2 librerías:

   
   

   1. Si vais con R: la scorecardModelUtils que está en CRAN

   2. Si vais con Python: está la opción (que también sirve para construir modelos) de OptBinning

Cuando se empieza un proyecto de este tipo, al poco tiempo, la empresa usuaria de los modelos se da cuenta del trabajo que realmente generan los modelos pero también de los beneficios uno de ellos es que muchas veces la historia no es como la de la gráfica anterior sino que es más bien del siguiente tipo (nótese que se mantienen las mismas escalas que en el gráfico anterior):

Y entonces cuando se observan estos datos y se comparan con los que nos ofreció el consultor, es cuando se toma conciencia de que algo no me han contado bien porque no quieren, bien porque no pueden o porque no sabían, ... El tema es que no se ve un claro ROI.

Y eso sólo a nivel modelo ¿Y a nivel variable? ¿Qué puede ocurrir? Pues, eso da para otro post y la cosa aterra aún más, sobre todo si usar alguna de esas variables extras tiene algún coste.

Finalmente indicar que en un sistema de monitoring esto no es más que la punta del iceberg, porque hay aún más cosas a monitorizar como ¿Nos están sirviendo las puntuaciones correctamente? ¿Las estamos generando adecuadamente? La cantidad de sorpresas que salen de esta parte cuando se mira un poco debajo de la alfombra da también para otro post, por tanto, que nadie crea que externalizar modelos ahorrará gente y dolores de cabeza, sino que, aunque en ocasiones puede ser necesario, hay que tener muy en cuenta, que va a ser sano que desde el momento uno se programe el monitoring, en caso contrario, el lío "cuando llegue yo", va ser monumental y solucionarlo, no va a ser muy agradable.

LO QUE NO TE CUENTA TU CONSULTOR PARTE III: El humo del Machine Learning de los primeros años

En un momento del tiempo llegó el Machine Learning para sustituir a la minería de datos a toda cosa anterior, algo a lo que está pasando con la IA, pero de un modo muy brusco se pasó de los análisis tradicionales a otros, en ocasiones algo más livianos, que sin ser malos del todo, hay que tener cierto cuidado.

A partir del 2014 se popularizan los svm, random forest y posteriormente irrumpirían los xgboost y demás modelos que me gusta llamarlos de analítica avanzada, junto con una pléyade de Data Scientist que básicamente metían vectores y sacaban ginis o roc altos y "a rular".

Así pues, el negocio de la consultoría se reinventa y ahora se pasó de estudiar la variable y de hablar de consistencia lógica de un modelo comprensible por el usuario, a algo que este no controlaba y lo peor de todo, el propio usuario lo prefería hasta el punto de llegar a escuchar afirmaciones del tipo "si lo entiendo es que es malo" y claro está, de nuevo surgen proyectos como setas tras unos días de lluvia en un bosque polaco.

Se generan modelos de esta nueva analítica y los más versados te presentaban gráficos del tipo importancia de las variables (normalizados o no como el siguiente):

Y los clientes, muchos sin entender pues ven que hay como una variable muy importante, otra sólo un poco importante, que cabría hasta la posibilidad de desecharla, pero nada como quién oye llover en muchas ocasiones porque le cuentan la historia de que si modelos débiles se unen para crear otro más fuete y si patatín o si patatán ...

Lo anterior se puede conseguir ejecutando por ejemplo el siguiente código:

library(dplyr)

library(scorecard)

library(ISLR)

library(ranger)

library(ggplot2)

data("Default")

df <- Default

df$default <- ifelse(df$default == "Yes", 1, 0)

df_list <- split_df(df, y = "default", ratios = c(0.6, 0.4), seed = 30)

label_list <- lapply(df_list, function(x) x$default)

df_train <- df_list$train

df_test <- df_list$test

# Modelo Random Forest con ranger

set.seed(123)

rf_model <- ranger(

formula = factor(default) ~ .,

data = df_train, importance = "impurity", probability = TRUE

)

imp <- data.frame(

Variable = names(rf_model$variable.importance),

Importancia = rf_model$variable.importance

)

ggplot(imp, aes(x = reorder(Variable, Importancia), y = Importancia)) +

geom_col(fill = "#0072B2") +

coord_flip() +

labs(

title = "Importancia de variables (Random Forest - ranger)",

x = "Variable",

y = "Importancia"

) +

theme_minimal()

Entonces ahora cabe hacerse la pregunta ¿Qué significa lo anterior? Lo anterior no es más que una medida de lo que se denomina impureza o cómo de mezcladas están las clases en un nodo de un árbol (o conjunto de árboles entremezclados como es este caso). Por lo que se tiene una medida bastante relativa y condicionada al dataset y al modelo en concreto plantead sin mucha más base teórica detrás, se elude el tema del signo y para alguien que no sabe, hasta te compra el modelo, máxime si le muestras el siguiente gráfico de desempeño:

pred_rf <- predict(rf_model, data = df_test)$predictions[,2]

perf_eva(label = df_test$default, pred = pred_rf,

title = "Curva Gini - Random Forest (ranger)",

show_plot = "roc",

binomial_metric = "gini")

Es decir, tenemos un AUC del 92.78% en test y un modelo de analítica avanzada donde las variables no se entienden bien cómo intervienen ¡Qué más se puede pedir!

Entonces ¿Qué ocurre si se trabaja desde la estadística tradicional y se plantea una regresión logística? Pues se obtiene la siguiente tabla característica:

m1 <- glm(default ~ ., data = df_train, family = binomial)

summary(m1)

Donde llama la atención que desde este punto de vista una variable como es la de income no es significativa y convendría eliminarse, además es contraintuitiva ¿Quiere decir que a mayor ingreso mayor probabilidad de impago? Lógicamente sale del modelo y quedaría tal como sigue:

Ahora parece que gana fuerza, desde un punto de vista del p-valor la variable student pero desde el sentido común y desde el siguiente gráfico asociado (visto en un post anterior), se tiene que univariantemente, a poco que se mire la variable el resultado es el siguiente:

Es decir, bajo el gráfico anterior el signo debería ser positivo, sin embargo resulta negativo, lo que va en contra del análisis univariante y es que como se puede observar a continuación, se observa que los que son estudiantes tienen un mayor nivel de balance y por tanto existe un grado de correlación hasta cierto punto peligrosa que conviene eliminar del modelo:

df_train %>%

group_by(student) %>%

summarise(m_balance = mean(balance))

Por tanto cabe plantear el modelo más simple y con sentido siguiente:

Y claro, si tú ibas con esto en el año 2018, te salían múltiples "expertiños" tanto de los consultores como del propio cliente diciendo que "si era una caca de la vaca", que si "no pagaban para cosas sencillas", que si ... Te habías currado el estudio de las variables, pero como ofreces algo que se entiende frente a un random forest dudoso, pues mejor tirarlo a la basura, prácticamente "te corrian a gorrazos", pues bien chavales ¿Cómo sería el performance de un sencillo modelo como el anterior? Si lo aplicamos sobre el test se obtiene lo siguiente:

Es decir, un AUC del 95.29, unos 3 puntos por encima y con sólo una variable, frente a 3 que usa un modelo de analítica avanzada como es el random forest, además se han descartado las variables por razones claras de sentido común y de negocio mediante su visualización.

Finalmente, y a pesar de este ejemplo, he visto en muchísimas ocasiones lo contrario, es decir que los modelos analítica avanzada suelen superar en la mayoría de los casos a los tradicionales en cuanto performance se refiere (si se manda al cuerno la interpretabilidad, claro está) y esta victoria se consigue cuando se introducen tablas en estos algoritmos y se trabaja lo mínimo las variables, pero cuando se hace algo que se quiere que tenga sentido, la diferencia ya no resulta tan aplastante, por tanto desde mi punto de vista, la irrupción en R y Python de estos modelos avanzados vienen a ser una caja de herramientas que debemos conocer los que nos dedicamos a esto, pero desde luego no nos dejemos impresionar por el humo de estas entelequias, pidamos contrastes con modelos más sencillos, pidamos análisis marginales de las contribuciones de las variables al modelo, veamos los EDAs que son fundamentales en cualquier proyecto y con toda una documentación clara y formal, entonces yo y cualquiera va aceptaremos un modelo de analítica avanzada, mientras tanto frente a los complicadores me aparto hacia la rama de los simplificadores y me complico eso sí, cuando sea necesario, no por modas.