1 Introducción

Nos encontramos ante un problema de aprendizaje supervisado, concretamente un problema de clasificación, pues los datos están etiquetados con clases, y las clases son en este caso son de tipo factor. En el presente trabajo se pretende construir modelos predictivos capaces de, recibiendo información sobre expresión génica y variables diversas sobre el paciente, predecir su estadio de braak, métrica empleada en la evaluación de pacientes con enfermedad de Alzheimer o Parkinson.

Ilustración del daño tisular y síntomas cognitivos en las distintas etapas braak, sacada de Petersen V., Mikkel (2017).



Primero vamos a cargar las librerias de turno.

library(parallel)
library(doParallel)
library(caret)
library(Boruta)
library(partykit)
library(pROC)
library(MLmetrics)
library(partykit)

  • Con las librerías parallel y doParallel podemos paralelizar los procedimientos aqui realizados, lo que aliviará los tiempos de espera en la ejecucion de ciertos comandos, sobre todo en los comandos de la libreria caret.

  • La librería caret será empleada para la generación de los modelos predictivos.

  • La librería Boruta se empleará para la selección de variables.

  • El paquete partykit será empleado para mostrar los árboles de decisiones en formato gráfico.

  • La librería pROC se usará para el análisis por curvas ROC.

  • La librería MLmetrics se usará para testear la presencia de overfitting en nuestros modelos.



Procedemos a establecer el nº de hilos que usaremos durante esta sesión

cores <- detectCores()-4 # Dejamos 4 hilos para el SO
registerDoParallel(cores)



Establecemos el espacio de trabajo de la sesión y cargamos los datos

setwd("C:/Users/Adam_/Desktop/Cosas master/Master Bioinformatica/2o Cuatrimestre/Machine Learning/Trabajos Finales ML/Clasificadores/")
dataset <- readRDS("expr.data.cov")

# Alternativamente, podemos descargarlo desde mi repositorio publico :)
# dataset2 <- readRDS(gzcon(url("https://github.com/gloknar/Adam-repo/blob/master/covariables?raw=true")))

# Echamos un vistazo a la naturaleza de las variables del dataset. Nótese que solo imprimimos las 15 últimas variables por cuestiones de legibilidad, ya que el resto de variables del dataset son numéricas y la información aportada sería saturante.
tail(sapply(dataset, class), n= 15L)
## ENSG00000143772 ENSG00000136848 ENSG00000164889 ENSG00000135631 ENSG00000119760 
##       "numeric"       "numeric"       "numeric"       "numeric"       "numeric" 
## ENSG00000236609           batch          gender             pmi             age 
##       "numeric"        "factor"        "factor"       "numeric"       "numeric" 
##            race         braaksc           cogdx            educ         ceradsc 
##        "factor"       "numeric"        "factor"        "factor"        "factor"



2 Preprocesado

Primero vamos a comprobar la presencia de valores NA mediante el siguiente bloque de código, el cual escanea todo el dataset, e imprime un mensaje por consola si encuentra un valor faltante. En caso contrario (como es aquí), no mostrará ningún mensaje. Si hubiese valores NA, se procedería a impugnar el dataset.

if(TRUE %in% apply(dataset, 2, is.na.data.frame)) {
  which(is.na(dataset) , arr.ind = T)
}



Ahora veremos si las variables numericas estan normalizadas y/o escaladas. Nótese que en este paso no tendremos en cuenta la columna nº70, pues es la clase braaksc, la cual es un factor, pero al estar codificada en números, R la identifica como una variable numérica.

variables_numericas <- as.numeric(which(sapply(dataset,class) == "numeric"))
variables_numericas <- variables_numericas[variables_numericas!=70]



Se observa que si bien los valores de expresión de los genes están centrados, no están escalados. Por otro lado, las variables numericas pmi y age (2 primeras cajas desde la derecha) no están centradas ni escaladas. Vamos a centrar y escalar dichas variables para mantener la media en 0 y a su vez evitar que algunas variables tengan más peso que otras en nuestros modelos predictivos por el mero hecho de presentar valores absolutos más grandes que el resto.

Un ejemplo de los sesgos que evitamos con el centrado y escalado de variables numéricas son los genes que codifican para factores de transcripción, pues según una comunicación personal del Dr. Javier Corral De la calle:

“Dichos genes suelen presentar niveles de expresión muy bajos en comparación con el resto, si bien tienen un efecto muy acusado en la regulación transcripcional, y por tanto, en el metabolismo de la célula.”


escalado_y_centrado <- preProcess(dataset[,variables_numericas], method = c("center", "scale"))
dataset_escalado <- predict(escalado_y_centrado,dataset[,variables_numericas])



A cotinuación eliminamos las variables numéricas altamente correlacionadas con el metodo findCorrelation del paquete caret. Para ello necesitamos crear una matriz de correlaciones que contenga exclusivamente las variables predictoras numéricas del dataset.

R <- cor(dataset_escalado)
caret::findCorrelation(
  R,
  cutoff = 0.75,
  verbose = TRUE,
  names = TRUE)
## Compare row 17  and column  13 with corr  0.754 
##   Means:  0.408 vs 0.31 so flagging column 17 
## All correlations <= 0.75
## [1] "ENSG00000105655"

En base a los resultados obtenidos, marcamos la columna 17 ENSG00000105655 para su eliminado, por estar altamente correlacionada con la columna 13 ENSG00000042062.

dataset_escalado <- dataset_escalado[,-17]


Seguidamente comprobamos en las columnas numéricas si hay variables con una varianza cercana a cero, y en caso afirmativo, eliminamos dichas columnas con el comando nzv.

nzv(dataset_escalado)
## integer(0)

Devuelve 0, por lo que no se encontraron columnas con varianza cercana a 0.



Nota antes de continuar con la selección de variables: si ploteamos el nº de individuos en función de la clase a la que pertenecen, se observa una acusada escasez de individuos que tienen diagnosticados los estadios braak 1 y 7, y esto puede dar problemas a la hora de entrenar los modelos, pues hay muy pocos individuos que representen a dichas clases, por lo que se procede a añadir los individuos con estadios braak 1 en el estadio braak 2, e ídem para los individuos con estadio braak 7, que pasan a formar parte del estadio braak 6


# Para individuos de la clase braak stage 1:
for (individuo in which(dataset$braaksc==1)) {dataset[individuo,"braaksc"] <- 2}

# Para individuos de la clase braak stage 7:
for (individuo in which(dataset$braaksc==7)) {dataset[individuo,"braaksc"] <- 6}



Una vez realizado este rebalanceo de clases, y teniendo en cuenta que la clase braaksc era reconocida (erróneamente) como una variable numérica, vamos a recodificarla como ordered factor con los siguientes comandos y comprobar que el cambio ha sido exitoso.

# Antes del type casting
tail(dataset$braaksc, n=0)
## numeric(0)

dataset$braaksc <- factor(dataset$braaksc,ordered = TRUE)

levels(dataset$braaksc) <-  c("estadios_1_2", "estadio_3", "estadio_4", "estadio_5","estadios_6_7")


# Después del type casting a factor ordenado
tail(dataset$braaksc, n=0)
## ordered(0)
## 5 Levels: estadios_1_2 < estadio_3 < estadio_4 < ... < estadios_6_7


Así mismo, los atributos cogdx y ceradsc son también ordered factors, a pesar de que son reconocidos como factors, por lo que vamos a castearlos a sus naturaleza adecuada.

Hay que tener cuidado con el atributo ceradsc, ya que en el caso de braaksc y cogdx, valores mayores de factor indican una clínica más severa; mientras que ceradsc está codificado la revés, de manera que valores menores del factor indican una clínica más severa. Para más información sobre la codificación de estos atributos, visite el RADC Research Resource Sharing Hub.

# Antes del type casting de la variable dependiente "cogdx" a factor ordenado
tail(dataset$cogdx, n = 0)
## factor(0)
## Levels: 1 2 3 4 5 6
dataset$cogdx <- factor(dataset$cogdx,ordered = TRUE)


# Después del type casting
tail(dataset$cogdx, n=0)
## ordered(0)
## Levels: 1 < 2 < 3 < 4 < 5 < 6
# ídem para ceradsc, pero inviertiendo el orden de los niveles:
tail(dataset$ceradsc, n = 0)
## factor(0)
## Levels: 1 2 3 4
dataset$ceradsc <- factor(dataset$ceradsc,ordered = TRUE)
levels(dataset$ceradsc) <-  c(4, 3, 2, 1)

# Después del type casting a factor ordenado
tail(dataset$ceradsc, n=0)
## ordered(0)
## Levels: 4 < 3 < 2 < 1


Finalmente, procedemos a añadir al dataset los factores que habíamos eliminado en el mismo durante su preprocesado.

# Obtenemos un vector con el nº de columna de los factores en el dataset original
columnas_factores <- as.numeric (which(sapply(dataset,class)== "factor"))
columnas_factores <- c(columnas_factores, which(colnames(dataset) %in% c("braaksc", "ceradsc", "cogdx"))) # añadimos factores ordenados

dataset_escalado <- cbind(dataset_escalado,dataset[,columnas_factores])

# Comprobamos que las columnas están bien formateadas. Por legibilidad, solo mostramos las 12 últimas.
tail(sapply(dataset_escalado, class), n = 12L)
## $ENSG00000135631
## [1] "numeric"
## 
## $ENSG00000119760
## [1] "numeric"
## 
## $ENSG00000236609
## [1] "numeric"
## 
## $pmi
## [1] "numeric"
## 
## $age
## [1] "numeric"
## 
## $batch
## [1] "factor"
## 
## $gender
## [1] "factor"
## 
## $race
## [1] "factor"
## 
## $educ
## [1] "factor"
## 
## $braaksc
## [1] "ordered" "factor" 
## 
## $cogdx
## [1] "ordered" "factor" 
## 
## $ceradsc
## [1] "ordered" "factor"



3 Selección de variables

Procedemos a la selección de variables. Usaremos el paquete Boruta para ello, dada su capacidad de mostrar visualmente la importancia de cada variable del dataset.

Boruta realiza una búsqueda de tipo Backwards mediante un algoritmo wrapper de Random Forest, y por defecto, usa un p-valor del 0.01 para confirmar cuándo una variable es estadísticamente importante o no.

# Establecemos la semilla de aleatoriedad para garantizar la reproducibilidad de
# esta sesion
set.seed(1)

# Ahora realizamos la busqueda de las variables importantes
boruta_output <- Boruta(braaksc ~ . , data=dataset_escalado, doTrace = 0)


La selección de variables tardó:

## Time difference of 15.06883 secs



Y la representación visual de dicha selección podemos verla en el siguiente gráfico:



Acto seguido nos quedamos con las variables significativamente importantes (i.e.“Confirmed”) y las imprimimos por consola

boruta_significativas <- names(boruta_output$finalDecision[boruta_output$finalDecision %in% "Confirmed"])  

print(boruta_significativas)   
## [1] "ENSG00000102547" "ENSG00000116337" "ENSG00000247982" "ENSG00000166012"
## [5] "ENSG00000149451" "ENSG00000135631" "age"             "cogdx"          
## [9] "ceradsc"



Posteriormente nos quedamos solo con las variables importantes y la clase, y las agrupamos en un nuevo dataframe llamado dataset_variables_seleccionadas.

variables_seleccionadas_boruta <- which(names(dataset_escalado) %in% boruta_significativas)
clase <- which(names(dataset_escalado) == "braaksc")

dataset_variables_seleccionadas <- dataset_escalado[,c(variables_seleccionadas_boruta,clase)]

sapply(dataset_variables_seleccionadas, class)
## $ENSG00000102547
## [1] "numeric"
## 
## $ENSG00000116337
## [1] "numeric"
## 
## $ENSG00000247982
## [1] "numeric"
## 
## $ENSG00000166012
## [1] "numeric"
## 
## $ENSG00000149451
## [1] "numeric"
## 
## $ENSG00000135631
## [1] "numeric"
## 
## $age
## [1] "numeric"
## 
## $cogdx
## [1] "ordered" "factor" 
## 
## $ceradsc
## [1] "ordered" "factor" 
## 
## $braaksc
## [1] "ordered" "factor"



4 Modelos predictivos

Vamos a hacer 3 modelos: uno hecho con J48, otro con una Red Neuronal (Perceptrón Multicapa) y otro con Regresión Logistica (modelo probabilístico). Pero antes, debemos partir el dataset. Lo haremos en un 80% para el training_set y un 20% para el test_set.

set.seed(1)
Indice_individuos_entrenamiento <- createDataPartition(dataset_variables_seleccionadas$braaksc,
                                  p = .80,
                                  list = FALSE,
                                  times = 1)

# Comprobamos que ha realizado la división deseada. Como en el comando anterior
# elegimos el parametro p = 0.80, estamos cogiendo el 80% de los individuos como
# set de entrenamiento. Nos debería dar en torno a 508 (635*0.8 = 508):
nrow(Indice_individuos_entrenamiento)
## [1] 511
# Creamos y comprobamos el set de entrenamiento:
training_set <- dataset_variables_seleccionadas[Indice_individuos_entrenamiento,]

cat("El set de entrenamiento tiene", nrow(training_set), "individuos y", ncol(training_set), "variables.")
## El set de entrenamiento tiene 511 individuos y 10 variables.
# Creamos y comprobamos el set de testeo:
test_set <- dataset_variables_seleccionadas[-Indice_individuos_entrenamiento,]

cat("El set de testeo tiene", nrow(test_set), "individuos y", ncol(test_set), "variables.")
## El set de testeo tiene 124 individuos y 10 variables.



A continuación establecemos las opciones de testeo que caret usará para la construcción de los clasificadores J48, MLP y RegLog (RegLog requiere un tratamiento ligeramente distinto, se verá en su correspondiente apartado). Para ellos crearemos primero un listado de semillas de aleatoriedad que usaremos más adelante para la creación de cada fold y la evaluación del modelo final:

set.seed(1)


# Iniciamos una lista vacía que almacenará las semillas de aleatoriedad con las
# que trabajaremos. Haremos 5 folds * 3 repeticiones = lista con 15 bolsillos +
# 1 para el modelo final (Lo pone en la documentación de ?trainControl() ).
# Debemos llenar todos los bolsillos con vectores con 15 numeros cada uno, menos
# en el último bolsillo de la lista semillas, que debe ser un único número.

semillas <- vector(mode = "list", length = 15)  


# Llenamos la lista/vector con semillas al azar, elegidas entre 1 y 1000.
# Creamos, por tanto, una lista de listas. Todos los bolsillos de la lista,
# menos el último, tendrán vectores con 16 números de longitud 
#
# (5 folds * 3 repeticiones) + 1 evaluación final = 16

semillas <- lapply(semillas, function(x) sample.int(1000, size= 300))

semillas[[16]] <- sample.int(1000, size= 1)

# Hemos usado lapply porque queremos que nos devuelva una lista/vector

str(semillas)
## List of 16
##  $ : int [1:300] 836 679 129 930 509 471 299 270 978 187 ...
##  $ : int [1:300] 270 751 209 338 609 736 584 565 568 743 ...
##  $ : int [1:300] 11 715 310 725 85 81 256 777 54 146 ...
##  $ : int [1:300] 201 487 211 5 465 329 694 492 182 89 ...
##  $ : int [1:300] 683 976 942 411 617 1 580 551 761 822 ...
##  $ : int [1:300] 165 744 966 759 953 248 690 117 435 751 ...
##  $ : int [1:300] 656 81 472 504 53 478 309 941 554 426 ...
##  $ : int [1:300] 442 225 344 519 307 815 114 429 723 126 ...
##  $ : int [1:300] 284 154 543 300 873 337 944 170 972 932 ...
##  $ : int [1:300] 608 742 464 811 568 601 929 367 349 810 ...
##  $ : int [1:300] 63 626 279 773 509 534 553 822 753 437 ...
##  $ : int [1:300] 370 68 589 503 934 26 166 18 312 765 ...
##  $ : int [1:300] 544 676 538 370 561 901 926 208 381 968 ...
##  $ : int [1:300] 769 560 391 32 551 398 139 534 90 861 ...
##  $ : int [1:300] 814 614 863 570 522 617 346 467 169 954 ...
##  $ : int 906


Con las semillas de aleatoriedad listas, establecemos el control del entrenamiento para nuestros modelos:

control_entrenamiento <- trainControl(method = "repeatedcv",
                           number=5,
                           repeats = 3,
                           seeds = semillas, #Poner "semillas_J48_MLP" cuando acabe
                           returnResamp = "final",
                           allowParallel = TRUE,
                           classProbs = TRUE,
                           summaryFunction = multiClassSummary, 
                           verboseIter = T,
                           search = "grid",
                           p = .8)

# Nota: La opción "summaryFunction = defaultSummary" solo permite obtener las
# métricas Accuracy y Kappa. Por otro lado,  "summaryFunction = multiClassSummary" 
# permite obtener el AUC (además del índice Kappa y la precisión) en problemas multiclase.


4.1 J48

4.1.1 Modelo con 8 predictores

Entrenamos el modelo con las variables seleccionadas por Boruta (8 variables predictoras + la clase) e intentando maximizar la métrica AUC. Además, por ser un problema multi-clase, prestaremos especia atención a la métrica Kappa de nuestros modelos.

El índice Kappa de Cohen mide la tasa de acuerdo entre las distintas evaluaciones realizadas por el clasificador de interés durante su entrenamiento, siendo un valor de 1 una tasa de acuerdo del 100%, y 0 una tasa del 0%. La figura a continuación detalla su interpretación:

Guía de interpretación índice kappa


Así mismo, el AUC o Area Under the Curve mide los distintos puntos de corte que podemos usar en nuestro modelo para clasificar las diversas clases del problema, y mide por tanto la tasa de verdaderos positivos vs falsos positivos detectados por el clasificador. La siguiente imagen trata de apoyar la explicación aquí aportada:


Nótese que el comando train nos permite seleccionar la métrica que queremos maximizar (o minimizar). Para todos nuestros clasificadores, maximizaremos la métrica AUC, ya que la usaremos posteriormente para evaluar el rendimiento de nuestros clasificadores.

set.seed(1)

modelo_J48 <- train(braaksc~.,data=training_set,
                  method="J48",
                  metric = "AUC",     # Maximize = TRUE by default
                  trControl = control_entrenamiento,
                  tuneLength = 10)
## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.
## Aggregating results
## Selecting tuning parameters
## Fitting C = 0.5, M = 10 on full training set


Los mejores hiperparámetros encontrados fueron:

modelo_J48$bestTune
# Cuadrícula de búsqueda de hiperparámetros óptimos de J48 (8 predictores)
ggplot(modelo_J48) + theme(legend.position = "top")


Ahora construimos el modelo final con los hiperparámetros óptimos (C= 0.5 y M = 10):

set.seed(1)

J48Grid <-  expand.grid(C = 0.5, M = 10)

modelo_J48 <- train(braaksc~.,data=training_set,
                  method= "J48",
                  metric = "AUC",     # Maximize = TRUE by default
                  trControl = control_entrenamiento,
                  tuneGrid = J48Grid)
## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.
## Aggregating results
## Fitting final model on full training set
# Se puede obtener información detallada sobre el modelo con los siguientes comandos:
# modelo_J48$results
# modelo_J48
cat("El área media bajo la curva ROC de este modelo fue:", mean(modelo_J48$resample$AUC), "\n")
## El área media bajo la curva ROC de este modelo fue: 0.7059865
cat("La precisión media de este modelo fue:", mean(modelo_J48$resample$Accuracy), "% \n")
## La precisión media de este modelo fue: 0.4554014 %
cat("El índice Kappa medio de este modelo fue:", mean(modelo_J48$resample$Kappa), "\n")
## El índice Kappa medio de este modelo fue: 0.2650073



Árbol generado:

(Aún si ampliásemos la imagen para verla a pantalla completa, podemos imaginar que un árbol de decisiones tan frondoso sería dificil de leer).


Visualizaremos mejor el árbol si lo mostramos con gráficos simples, tal que así:

## J48 pruned tree
## ------------------
## 
## ceradsc.L <= -0.223607
## |   ceradsc.L <= -0.67082
## |   |   cogdx.L <= -0.119523
## |   |   |   age <= -0.223271: estadio_5 (14.0/6.0)
## |   |   |   age > -0.223271
## |   |   |   |   cogdx.L <= -0.597614: estadios_6_7 (13.0/7.0)
## |   |   |   |   cogdx.L > -0.597614
## |   |   |   |   |   ENSG00000102547 <= -0.153971: estadio_5 (11.0/4.0)
## |   |   |   |   |   ENSG00000102547 > -0.153971: estadios_6_7 (10.0/4.0)
## |   |   cogdx.L > -0.119523
## |   |   |   ENSG00000166012 <= -0.505501: estadios_6_7 (52.0/8.0)
## |   |   |   ENSG00000166012 > -0.505501
## |   |   |   |   ENSG00000149451 <= 0.010508: estadios_6_7 (28.0/13.0)
## |   |   |   |   ENSG00000149451 > 0.010508
## |   |   |   |   |   ENSG00000247982 <= 0.062957: estadios_6_7 (11.0/5.0)
## |   |   |   |   |   ENSG00000247982 > 0.062957: estadio_5 (10.0/3.0)
## |   ceradsc.L > -0.67082
## |   |   age <= -0.752769
## |   |   |   ENSG00000247982 <= -0.245729: estadio_3 (16.0/9.0)
## |   |   |   ENSG00000247982 > -0.245729: estadio_5 (11.0/5.0)
## |   |   age > -0.752769
## |   |   |   cogdx.Q <= -0.109109
## |   |   |   |   ENSG00000116337 <= -0.948559: estadio_5 (20.0/3.0)
## |   |   |   |   ENSG00000116337 > -0.948559
## |   |   |   |   |   ENSG00000102547 <= -0.654521: estadio_5 (16.0/3.0)
## |   |   |   |   |   ENSG00000102547 > -0.654521
## |   |   |   |   |   |   cogdx^5 <= 0.31497
## |   |   |   |   |   |   |   ENSG00000135631 <= 0.04203: estadio_5 (13.0/6.0)
## |   |   |   |   |   |   |   ENSG00000135631 > 0.04203: estadio_4 (16.0/8.0)
## |   |   |   |   |   |   cogdx^5 > 0.31497
## |   |   |   |   |   |   |   ENSG00000135631 <= 0.04203: estadio_4 (14.0/7.0)
## |   |   |   |   |   |   |   ENSG00000135631 > 0.04203: estadios_6_7 (25.0/9.0)
## |   |   |   cogdx.Q > -0.109109
## |   |   |   |   ENSG00000149451 <= -0.364543: estadio_4 (13.0/6.0)
## |   |   |   |   ENSG00000149451 > -0.364543: estadio_5 (32.0/12.0)
## ceradsc.L > -0.223607
## |   age <= 0.21828
## |   |   cogdx^4 <= -0.566947: estadio_4 (27.0/12.0)
## |   |   cogdx^4 > -0.566947
## |   |   |   ENSG00000149451 <= 0.502476
## |   |   |   |   cogdx.L <= -0.597614
## |   |   |   |   |   ENSG00000102547 <= 0.549232
## |   |   |   |   |   |   ENSG00000116337 <= -0.36724: estadios_1_2 (16.0/6.0)
## |   |   |   |   |   |   ENSG00000116337 > -0.36724: estadio_4 (17.0/10.0)
## |   |   |   |   |   ENSG00000102547 > 0.549232: estadio_4 (10.0/1.0)
## |   |   |   |   cogdx.L > -0.597614: estadios_1_2 (10.0/4.0)
## |   |   |   ENSG00000149451 > 0.502476: estadios_1_2 (16.0/5.0)
## |   age > 0.21828
## |   |   ENSG00000149451 <= -0.558676: estadio_4 (38.0/19.0)
## |   |   ENSG00000149451 > -0.558676
## |   |   |   cogdx.L <= -0.119523: estadio_4 (36.0/19.0)
## |   |   |   cogdx.L > -0.119523: estadio_5 (16.0/5.0)
## 
## Number of Leaves  :  27
## 
## Size of the tree :   53


Por último, comprobamos si hay overfitting comparando la precisión media del modelo en el set de entrenamiento y en el de test:

# Precisión en el set de entrenamiento
A <- MLmetrics::Accuracy(predict(modelo_J48, training_set[,-10]), training_set[[10]])

# Precisión en el set de test
B <- MLmetrics::Accuracy(predict(modelo_J48, test_set[,-10]), test_set[[10]])


print(paste0("La diferencia entre la precisión en el set de entrenamiento (", A, "%) y en la del set de test (", B, "%) es: ", A - B,"%" ))
## [1] "La diferencia entre la precisión en el set de entrenamiento (0.610567514677104%) y en la del set de test (0.387096774193548%) es: 0.223470740483555%"

Podemos observar que hay una diferencia importante de rendimiento entre el set de entrenamiento y el de test (~22%), por lo que podemos asegurar que se da overfitting en este modelo.



4.1.2 Modelo con 4 predictores

Puesto que con 8 predictores se genera un árbol bastante lioso, vamos a probar a volver a entrenar el modelo, pero esta vez con sólo 4 variables predictoras más la clase. Recordemos que las 4 variables más importantes son, en orden decreciente: ceradsc, cogdx, age y ENSG00000135631. Con esta simplificación del modelo esperamos mejorar su capacidad para generalizar, además de mejorar sustancialmente la comprensibilidad del mismo.

columnas_J48_2 <- vector()

for (nombre in c("ceradsc", "cogdx", "age", "ENSG00000135631", "braaksc")) {
  x <- which(colnames(training_set) == nombre)
  columnas_J48_2 <- c(columnas_J48_2, x)
}

# Comprobamos que se realizó correctamente el bucle
colnames(training_set)[columnas_J48_2]
## [1] "ceradsc"         "cogdx"           "age"             "ENSG00000135631"
## [5] "braaksc"


Construimos el modelo y buscamos los mejores hiperparámetros:

set.seed(1)
modelo_J48_2 <- train(braaksc~.,data=training_set[,columnas_J48_2],
                    method="J48",
                    metric = "AUC",      # Maximize = TRUE by default
                    tuneLength = 10,
                    trControl = control_entrenamiento)
## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.
## Aggregating results
## Selecting tuning parameters
## Fitting C = 0.337, M = 10 on full training set


La mejor configuración de hiperparámetros encontrada fue:

modelo_J48_2$bestTune
# Cuadrícula de búsqueda de hiperparámetros óptimos de J48 (4 predictores)
ggplot(modelo_J48_2) + theme(legend.position = "top")


Ahora construimos el modelo final con los hiperparámetros óptimos (C= 0.33 y M = 10):

set.seed(1)

J48_2Grid <-  expand.grid(C = 0.33, M = 10)

modelo_J48_2 <- train(braaksc~.,data=training_set[,columnas_J48_2],
                  method= "J48",
                  metric = "AUC",     # Maximize = TRUE by default
                  trControl = control_entrenamiento,
                  tuneGrid = J48_2Grid)
## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.
## Aggregating results
## Fitting final model on full training set
# Se puede obtener información detallada sobre el modelo con los siguientes comandos:
# modelo_J48_2$results
# modelo_J48_2
cat("El área media bajo la curva ROC de este modelo fue:", mean(modelo_J48_2$resample$AUC), "\n")
## El área media bajo la curva ROC de este modelo fue: 0.7230482
cat("La precisión media de este modelo fue:", mean(modelo_J48_2$resample$Accuracy), "% \n")
## La precisión media de este modelo fue: 0.4692049 %
cat("El índice Kappa medio de este modelo fue:", mean(modelo_J48_2$resample$Kappa), "\n")
## El índice Kappa medio de este modelo fue: 0.2851677



Árbol generado:

Si bien la precisión de este segundo modelo parece ser similar a la de J48 con 8 predictores (lo contrastaremos más adelante), también es cierto que este árbol es mucho más fácil de intepretar.


Visualizaremos mejor el árbol si lo mostramos con gráficos simples, tal que así:

## J48 pruned tree
## ------------------
## 
## ceradsc.L <= -0.223607
## |   ceradsc.L <= -0.67082
## |   |   cogdx.L <= -0.119523
## |   |   |   age <= -0.223271: estadio_5 (14.0/6.0)
## |   |   |   age > -0.223271
## |   |   |   |   ENSG00000135631 <= 0.031108: estadio_5 (19.0/8.0)
## |   |   |   |   ENSG00000135631 > 0.031108: estadios_6_7 (15.0/6.0)
## |   |   cogdx.L > -0.119523: estadios_6_7 (101.0/34.0)
## |   ceradsc.L > -0.67082
## |   |   age <= -0.752769
## |   |   |   ENSG00000135631 <= -0.518991: estadio_3 (12.0/7.0)
## |   |   |   ENSG00000135631 > -0.518991: estadio_5 (15.0/8.0)
## |   |   age > -0.752769
## |   |   |   ENSG00000135631 <= 0.054551: estadio_5 (75.0/28.0)
## |   |   |   ENSG00000135631 > 0.054551
## |   |   |   |   cogdx.L <= -0.119523
## |   |   |   |   |   ENSG00000135631 <= 0.637668: estadio_4 (24.0/11.0)
## |   |   |   |   |   ENSG00000135631 > 0.637668: estadio_5 (18.0/8.0)
## |   |   |   |   cogdx.L > -0.119523: estadios_6_7 (32.0/15.0)
## ceradsc.L > -0.223607
## |   age <= 0.21828
## |   |   cogdx^4 <= -0.566947: estadio_4 (27.0/12.0)
## |   |   cogdx^4 > -0.566947
## |   |   |   cogdx.L <= -0.597614
## |   |   |   |   ENSG00000135631 <= -0.799497: estadios_1_2 (18.0/5.0)
## |   |   |   |   ENSG00000135631 > -0.799497: estadio_4 (40.0/21.0)
## |   |   |   cogdx.L > -0.597614: estadios_1_2 (11.0/4.0)
## |   age > 0.21828
## |   |   cogdx.C <= -0.298142
## |   |   |   cogdx.L <= -0.358569
## |   |   |   |   ceradsc.L <= 0.223607: estadio_4 (12.0/6.0)
## |   |   |   |   ceradsc.L > 0.223607
## |   |   |   |   |   ENSG00000135631 <= -0.025175: estadio_3 (10.0/6.0)
## |   |   |   |   |   ENSG00000135631 > -0.025175: estadio_4 (11.0/5.0)
## |   |   |   cogdx.L > -0.358569
## |   |   |   |   ENSG00000135631 <= -0.016768: estadio_4 (14.0/7.0)
## |   |   |   |   ENSG00000135631 > -0.016768: estadio_5 (11.0/5.0)
## |   |   cogdx.C > -0.298142: estadio_4 (32.0/13.0)
## 
## Number of Leaves  :  20
## 
## Size of the tree :   39


Por último, comprobamos si hay overfitting comparando la precisión media del modelo en el set de entrenamiento y en el de test:

# Precisión en el set de entrenamiento
A <- MLmetrics::Accuracy(predict(modelo_J48_2, training_set[,c(6, 7, 8, 9)]), training_set[[10]])

# Precisión en el set de test
B <- MLmetrics::Accuracy(predict(modelo_J48_2, test_set[,c(6, 7, 8, 9)]), test_set[[10]])

print(paste0("La diferencia entre la precisión en el set de entrenamiento (", A, "%) y en la del set de test (", B, "%) es: ", A - B,"%" ))
## [1] "La diferencia entre la precisión en el set de entrenamiento (0.579256360078278%) y en la del set de test (0.395161290322581%) es: 0.184095069755697%"

La precisión en el set de entrenamiento es mejor que en la de test, pero no demasiado (~18% de diferencia), así que no podemos asegurar con total certeza que haya overfitting.



4.2 MLP

Entrenamos el modelo con las variables seleccionadas por Boruta (8 variables predictoras + la clase).

set.seed(1)
modelo_MLP <- train(braaksc~.,data=training_set,
                method="mlp",
                metric = "AUC",      # Maximize = TRUE by default
                tuneLength = 10,
                trControl = control_entrenamiento)
## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.
## Aggregating results
## Selecting tuning parameters
## Fitting size = 5 on full training set


La mejor configuración de hiperparámetros encontrada fue:

modelo_MLP$bestTune
# Cuadrícula de búsqueda de hiperparámetros óptimos de MLP
ggplot(modelo_MLP) + theme(legend.position = "top")


Ahora construimos el modelo final con los hiperparámetros óptimos (Hidden units = 3):

set.seed(1)

MLPGrid <-  expand.grid(size = 3)

modelo_MLP <- train(braaksc~.,data=training_set,
                method="mlp",
                metric = "AUC",      # Maximize = TRUE by default
                trControl = control_entrenamiento,
                tuneGrid = MLPGrid)
## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.
## Aggregating results
## Fitting final model on full training set
# Se puede obtener información detallada sobre el modelo con los siguientes comandos:
# modelo_MLP$results
# modelo_MLP
cat("El área media bajo la curva ROC de este modelo fue:", mean(modelo_MLP$resample$AUC), "\n")
## El área media bajo la curva ROC de este modelo fue: 0.7484993
cat("La precisión media de este modelo fue:", mean(modelo_MLP$resample$Accuracy), "% \n")
## La precisión media de este modelo fue: 0.4580562 %
cat("El índice Kappa medio de este modelo fue:", mean(modelo_MLP$resample$Kappa), "\n")
## El índice Kappa medio de este modelo fue: 0.2595578



Por último, comprobamos si hay overfitting comparando la precisión media del modelo en el set de entrenamiento y en el de test:

# Precisión en el set de entrenamiento
A <- MLmetrics::Accuracy(predict(modelo_MLP, training_set[,-10]), training_set[[10]])

# Precisión en el set de test
B <- MLmetrics::Accuracy(predict(modelo_MLP, test_set[,-10]), test_set[[10]])

print(paste0("La diferencia entre la precisión en el set de entrenamiento (", A, "%) y en la del set de test (", B, "%) es: ", A - B,"%" ))
## [1] "La diferencia entre la precisión en el set de entrenamiento (0.532289628180039%) y en la del set de test (0.419354838709677%) es: 0.112934789470362%"

La precisión en el set de entrenamiento es ligeramente mejor que en la de test (~11% de diferencia), por lo que puede que haya overfitting, pero no estamos seguros.



4.3 Regresión logística

Al igual que el primer modelo de J48 y el de MLP, esta regresión también se realizará con 8 predictores más la clase. Nótese que RegLog requiere de 300 semillas de aleatoriedad durante su entrenamiento (gradilla de búsqueda de 10 valores * 3 hiperparámetros).


set.seed(1)
modelo_RegLog <- train(braaksc~.,data=training_set,
                method="regLogistic",
                metric = "AUC",      # Maximize = TRUE by default
                tuneLength = 10,
                trControl = control_entrenamiento)
## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.
## Aggregating results
## Selecting tuning parameters
## Fitting cost = 0.5, loss = L1, epsilon = 0.01 on full training set


La mejor configuración de hiperparámetros encontrada fue:

modelo_RegLog$bestTune
# Cuadrícula de búsqueda de hiperparámetros óptimos de RegLog
ggplot(modelo_RegLog) + theme(legend.position = "top")


Ahora construimos el modelo final con los hiperparámetros óptimos (cost= 0.5, y loss = L1 y epsilon = 0.01):

set.seed(1)

RegLogGrid <-  expand.grid(cost = 0.5, loss = "L1", epsilon = 0.01)


modelo_RegLog <- train(braaksc~.,data=training_set,
                method="regLogistic",
                metric = "AUC",      # Maximize = TRUE by default
                trControl = control_entrenamiento,
                tuneGrid = RegLogGrid)
## Warning in nominalTrainWorkflow(x = x, y = y, wts = weights, info = trainInfo, :
## There were missing values in resampled performance measures.
## Aggregating results
## Fitting final model on full training set
# Se puede obtener información detallada sobre el modelo con los siguientes comandos:
# modelo_RegLog$results
# modelo_RegLog
cat("El área media bajo la curva ROC de este modelo fue:", mean(modelo_RegLog$resample$AUC), "\n")
## El área media bajo la curva ROC de este modelo fue: 0.7655036
cat("La precisión media de este modelo fue:", mean(modelo_RegLog$resample$Accuracy), "% \n")
## La precisión media de este modelo fue: 0.5017857 %
cat("El índice Kappa medio de este modelo fue:", mean(modelo_RegLog$resample$Kappa), "\n")
## El índice Kappa medio de este modelo fue: 0.3161772


Por último, comprobamos si hay overfitting comparando la precisión media del modelo en el set de entrenamiento y en el de test:

# Precisión en el set de entrenamiento
A <- MLmetrics::Accuracy(predict(modelo_RegLog, training_set[,-10]), training_set[[10]])

# Precisión en el set de test
B <- MLmetrics::Accuracy(predict(modelo_RegLog, test_set[,-10]), test_set[[10]])

print(paste0("La diferencia entre la precisión en el set de entrenamiento (", A, "%) y en la del set de test (", B, "%) es: ", A - B,"%" ))
## [1] "La diferencia entre la precisión en el set de entrenamiento (0.547945205479452%) y en la del set de test (0.435483870967742%) es: 0.11246133451171%"

La precisión en el set de entrenamiento es ligeramente mejor que en la de test (~11% de diferencia), por lo que puede que haya overfitting, pero no estamos seguros.



5 Análisis estadístico del rendimiento de los clasificadores

La Regresión Logística con 8 predictores parece ser el clasificador que mayor precisión, índice Kappa y AUC devuelve de los 4 modelos aquí presentes (téngase en cuenta que por ser un problema multi-clase, lo ideal es prestar especial atención al índice Kappa, en lugar de a la métrica de precisión o Accuracy).

Modelo Precisión Kappa AUC
J48 (8 predictores) 0.4554014% 0.2650073 0.7059865
J48 (4 predictores) 0.4692049% 0.2851677 0.7230482
MLP 0.4580562% 0.2595578 0.7484993
RegLog 0.5017857% 0.3161772 0.7655036

Nota: Obsérvese que tal como anticipamos, las métricas del modelo J48 mejoraron ligeramente al simplificarlo 😄


No obstante, antes de afirmar nada, procederemos a realizar un contraste de medias para la métrica Kappa de los modelos en busca de diferencias estadísticamente significativas. Puesto que durante el entrenamiento de los modelos, estos se evaluaron mediante el método repeatedcv con 5 pliegues y 3 repeticiones, cada modelo ha generado (5 pliegues*3 repeticiones + 1 evaluación del modelo final =) 16 repeticiones de las métricas. Con los siguientes comando de R podemos ver un resumen de las 15 evaluaciones realizadas durante el entrenamiento de los modelos, y también podemos examinarlas individualmente, pliegue a pliegue, repetición a repetición. Nota: Los valores en la tabla superior correspoden a la evaluación de los modelos finales.


modelos <- list(J48=modelo_J48, J48_2=modelo_J48_2, MLP=modelo_MLP, RegLog=modelo_RegLog)
datos_muestreo <- resamples(modelos)

# Como se calcularon muchas métricas, nos quedamos con las que nos interesan, que son la precisión, el Kappa y el AUC de nuestros 4 modelos
datos_muestreo$values  <- datos_muestreo$values[,c(1:4,16:18,30:32,44:46)]

# Adicionalmente, podemos analizar individualmente cada fold con este comando:
datos_muestreo$values


Antes de realizar el contraste de medias por test t de Student o ANOVA, primero comprobaremos si se cumplen en las muestras las condiciones de normalidad y homocedasticidad. Para la normalidad, usamos el test de normalidad de Saphiro-Wilk por tener un tamaño muestral n = 15 <=> n <= 50. En caso de no cumplirse alguna de las 2 condiciones, recurriremos a un test no paramétrico equivalente al t de Student.

metricas_kappa <- datos_muestreo$values[,c("J48~Kappa","J48_2~Kappa","MLP~Kappa","RegLog~Kappa")]


metricas_accuracy <- datos_muestreo$values[,c("J48~Accuracy","J48_2~Accuracy","MLP~Accuracy","RegLog~Accuracy")]


metricas_auc <- datos_muestreo$values[,c("J48~AUC","J48_2~AUC","MLP~AUC","RegLog~AUC")]


Procedemos al contraste de la normalidad por test de Shapiro-Wilk, cuya hipótesis nula H0 establece que los datos son normales:

for(i in 1:ncol(metricas_kappa)) {
  print(shapiro.test(metricas_kappa[,i]))
}
## 
##  Shapiro-Wilk normality test
## 
## data:  metricas_kappa[, i]
## W = 0.97419, p-value = 0.9145
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  metricas_kappa[, i]
## W = 0.89548, p-value = 0.0812
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  metricas_kappa[, i]
## W = 0.95146, p-value = 0.5478
## 
## 
##  Shapiro-Wilk normality test
## 
## data:  metricas_kappa[, i]
## W = 0.93767, p-value = 0.354

Si aplicamos la corrección de Bonferroni al alfa, tenemos un valor crítico de 0.05/4 = 0.0125.

Vemos que las métricas Kappa de los modelos son normales, ya que todos nuestros p-valores son mayores a 0.0125 y por tanto no podemos rechazar la hipótesis nula.


A continuación hacemos el contraste para la igualdad de varianzas de los índices Kappa. Como son 4 poblaciones normales e independientes, usaremos el test de esfericidad de Bartlett, cuya hipótesis nula H0 establece que se cumple la condición de homocedasticidad en las muestras.

bartlett.test(list(metricas_kappa$`J48~Kappa`, metricas_kappa$`J48_2~Kappa`, metricas_kappa$`MLP~Kappa`, metricas_kappa$`RegLog~Kappa`))
## 
##  Bartlett test of homogeneity of variances
## 
## data:  list(metricas_kappa$`J48~Kappa`, metricas_kappa$`J48_2~Kappa`, metricas_kappa$`MLP~Kappa`, metricas_kappa$`RegLog~Kappa`)
## Bartlett's K-squared = 1.245, df = 3, p-value = 0.7422

Puesto que nuestro p-valor es mayor a 0.05, no rechazamos la hipótesis nula. Vemos por tanto que las métricas Kappa de los modelos tienen varianzas iguales.


Ya que los índices Kappa de nuestros clasificadores son normales y tienen varianzas iguales, podemos aplicarles el test t de Student o un ANOVA. Vamos a mostrar primero cómo se haría con un ANOVA, y luego mostraremos cómo se haría con múltiples tests t de Student:

matriz_anova_kappa <- as.data.frame(metricas_kappa$`J48~Kappa`)
matriz_anova_kappa <- cbind(matriz_anova_kappa, rep("J48", 15))
colnames(matriz_anova_kappa) <- c("Kappa","Modelo")


variable_aux <- as.data.frame(metricas_kappa$`J48_2~Kappa`)
variable_aux <- cbind(variable_aux, rep("J48_2", 15))
colnames(variable_aux) <- c("Kappa","Modelo")


variable_aux2 <- as.data.frame(metricas_kappa$`MLP~Kappa`)
variable_aux2 <- cbind(variable_aux2, rep("MLP", 15))
colnames(variable_aux2) <- c("Kappa","Modelo")



variable_aux3 <- as.data.frame(metricas_kappa$`RegLog~Kappa`)
variable_aux3 <- cbind(variable_aux3, rep("RegLog", 15))
colnames(variable_aux3) <- c("Kappa","Modelo")



matriz_anova_kappa <- rbind2(matriz_anova_kappa, variable_aux)
matriz_anova_kappa <- rbind2(matriz_anova_kappa, variable_aux2)
matriz_anova_kappa <- rbind2(matriz_anova_kappa, variable_aux3)
matriz_anova_kappa$Modelo <- as.factor(matriz_anova_kappa$Modelo)

str(matriz_anova_kappa)
## 'data.frame':    60 obs. of  2 variables:
##  $ Kappa : num  0.226 0.202 0.258 0.301 0.291 ...
##  $ Modelo: Factor w/ 4 levels "J48","J48_2",..: 1 1 1 1 1 1 1 1 1 1 ...
anova_kappa <- aov(Kappa~Modelo, data = matriz_anova_kappa)
summary.aov(anova_kappa)
##             Df  Sum Sq  Mean Sq F value Pr(>F)  
## Modelo       3 0.02954 0.009847   3.179 0.0309 *
## Residuals   56 0.17345 0.003097                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


El ANOVA detectó diferencias significativas entre los índice Kappa de los diversos clasificadores. Procedemos con el test de Tukey para ahondar en los resultados del ANOVA:

TukeyHSD(anova_kappa)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Kappa ~ Modelo, data = matriz_anova_kappa)
## 
## $Modelo
##                      diff          lwr        upr     p adj
## J48_2-J48     0.020160395 -0.033648931 0.07396972 0.7545222
## MLP-J48      -0.005449476 -0.059258802 0.04835985 0.9931774
## RegLog-J48    0.051169930 -0.002639396 0.10497926 0.0679635
## MLP-J48_2    -0.025609871 -0.079419197 0.02819946 0.5917264
## RegLog-J48_2  0.031009535 -0.022799791 0.08481886 0.4291322
## RegLog-MLP    0.056619406  0.002810079 0.11042873 0.0355456

Hay diferencias estadísticamente significativas entre el rendimiento a nivel de índice Kappa de los clasificadores MLP y RegLog. En cuanto a la comparación RegLog vs J48 (8 predictores), presenta un p-valor ajustado cercano al punto crítico, pero superior al mismo, por lo que aceptamos parcialmente la hipótesis nula (es posible que haya diferencias significativas entre dicha pareja de modelos, pero nuestros resultados actuales carecen de la resolución necesaria para afirmarlo con seguridad).




A continuación mostramos cómo se realiza el contraste de medias con múltiples tests t de Student:

t.test(metricas_kappa$`J48~Kappa`, metricas_kappa$`J48_2~Kappa`, alternative = "two.sided", var.equal = TRUE) # J48 (8 pred.) vs J48 (4 pred.)
## 
##  Two Sample t-test
## 
## data:  metricas_kappa$`J48~Kappa` and metricas_kappa$`J48_2~Kappa`
## t = -0.96581, df = 28, p-value = 0.3424
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.06291882  0.02259803
## sample estimates:
## mean of x mean of y 
## 0.2650073 0.2851677
t.test(metricas_kappa$`J48~Kappa`, metricas_kappa$`MLP~Kappa`, alternative = "two.sided", var.equal = TRUE) # J48 (8 pred.) vs MLP
## 
##  Two Sample t-test
## 
## data:  metricas_kappa$`J48~Kappa` and metricas_kappa$`MLP~Kappa`
## t = 0.28953, df = 28, p-value = 0.7743
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.03310544  0.04400439
## sample estimates:
## mean of x mean of y 
## 0.2650073 0.2595578
t.test(metricas_kappa$`J48~Kappa`, metricas_kappa$`RegLog~Kappa`, alternative = "two.sided", var.equal = TRUE) # J48 (8 pred.) vs RegLog
## 
##  Two Sample t-test
## 
## data:  metricas_kappa$`J48~Kappa` and metricas_kappa$`RegLog~Kappa`
## t = -2.7462, df = 28, p-value = 0.01042
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.08933724 -0.01300262
## sample estimates:
## mean of x mean of y 
## 0.2650073 0.3161772
t.test(metricas_kappa$`J48_2~Kappa`, metricas_kappa$`MLP~Kappa`, alternative = "two.sided", var.equal = TRUE) # J48 (4 pred.) vs MLP
## 
##  Two Sample t-test
## 
## data:  metricas_kappa$`J48_2~Kappa` and metricas_kappa$`MLP~Kappa`
## t = 1.1704, df = 28, p-value = 0.2517
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.01921042  0.07043016
## sample estimates:
## mean of x mean of y 
## 0.2851677 0.2595578
t.test(metricas_kappa$`J48_2~Kappa`, metricas_kappa$`RegLog~Kappa`, alternative = "two.sided", var.equal = TRUE) # J48 (4 pred.) vs RegLog
## 
##  Two Sample t-test
## 
## data:  metricas_kappa$`J48_2~Kappa` and metricas_kappa$`RegLog~Kappa`
## t = -1.4278, df = 28, p-value = 0.1644
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.07549684  0.01347777
## sample estimates:
## mean of x mean of y 
## 0.2851677 0.3161772
t.test(metricas_kappa$`MLP~Kappa`, metricas_kappa$`RegLog~Kappa`, alternative = "two.sided", var.equal = TRUE) # MLP vs RegLog
## 
##  Two Sample t-test
## 
## data:  metricas_kappa$`MLP~Kappa` and metricas_kappa$`RegLog~Kappa`
## t = -2.8663, df = 28, p-value = 0.007799
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.09708320 -0.01615561
## sample estimates:
## mean of x mean of y 
## 0.2595578 0.3161772

Como hemos realizado múltiples tests, vamos a aplicarles a sus p-valores la corrección de Benjamini-Hochberg ( Yoav Benjamini & Yosef Hochberg, 1995 ). Se eligió esta corrección y no la de Bonferroni por ser más potente y fiable. Este procedimiento consiste en ordenar los p-valores de los tests de menor a mayor, y rechazar la hipótesis nula de todos los tests cuyo p-valor cumplan la siguiente condición:

\[\text{p-valor}_{test} < \frac{\text{rango del test}}{\text{nº total de tests}} *0.05\]

Nota: el 0.05 de la fórmula surge de elegir el FDR al 5%.


Resultados de aplicar la corrección de Benjamini-Hochberg.


De estos resultados podemos inferir las siguientes conclusiones:

  • La Regresión Logística ofrece un rendimiento significativamente distinto (=superior) al de MLP y posiblemente al de J48 con 8 predictores también (pero este test presenta un p-valor cercano al punto de corte, por lo que rechazamos parcialmente la hipótesis nula).
  • No existen diferencias estadísticamente significativas entre el índice Kappa del resto de modelos.



5.1 Conocimiento adicional inferido

Antes de avanzar al siguiente punto, pongamos en perspectiva la precisión de nuestros modelos. Según el Dr. Ian Witten en su MOOC Data mining with Weka, la precisión mínima posible para cualquier problema de clasificación se puede calcular como:

\[Accuracy_{min}=\frac{\text{major class}}{n}\]

siendo major class el nº de individuos de la clase más frecuente y n el nº de individuos en el dataset (esta es la base del algoritmo ZeroR, por cierto).

En este caso, la precisión base del problema es \(\frac{209}{635}= 0,3291338\) ó \(33\%\) (recordemos que la clase mayoritaria es estadio braak 5, con 209 individuos). Si observamos los boxplots de la precisión de los clasificadores, podemos ver que todos ellos son significativamente mejores que dicho valor umbral. No obstante, como se explico anteriormente, el foco de este estudio es el índice Kappa de los distintos modelos, procedimiento el cual ya se ha realizado.


Adicionalmente, con el comando varImp() podemos observar la relevancia que le dan los modelos a los predictores, para cada clase. Se observa que entre los diversos modelos hay consenso en:

  • A la hora de predecir las clases 1 y 2, la variable más importante fue la edad.
  • Para la clase 3, estas fueron ceradsc y la edad.
  • Para la clase 4, fue ceradsc y (en menor medida) cogdx y la edad.
  • Para la clase 5, el predictor más relevante fue ceradsc.
  • Para las clases 6 y 7, estas fueron ceradsc y la edad.






6 Evaluación de los modelos por métrica AUC

Esta evaluación se puede hacer de manera similar a como la hemos realizado con la métrica Kappa, o bien usar el comando multiclass.roc del paquete pROC. A continuación se demuestra el procedimiento con el susodicho comando:

6.1 J48

prediccion_J48 <- as.numeric(predict(modelo_J48, training_set, type = 'raw'))
AUC_J48 <- multiclass.roc(training_set$braaksc, prediccion_J48)
AUC_J48$auc
## Multi-class area under the curve: 0.8198

Para el modelo de J48 con 8 predictores, devuelve un AUC de 0.8367.


prediccion_J48_2 <- as.numeric(predict(modelo_J48_2, training_set, type = 'raw'))

AUC_J48_2 <- multiclass.roc(training_set$braaksc, prediccion_J48_2)
AUC_J48_2$auc
## Multi-class area under the curve: 0.7997

Para el modelo de J48 con 4 predictores, devuelve un AUC de 0.7549.



6.2 MLP

prediccion_MLP <- as.numeric(predict(modelo_MLP, training_set, type = 'raw'))

AUC_MLP <- multiclass.roc(training_set$braaksc, prediccion_MLP)
AUC_MLP$auc
## Multi-class area under the curve: 0.7941

Para el modelo de MLP con 8 predictores, devuelve un AUC de 0.7998.



6.3 Regresión Logística

set.seed(1)
prediccion_RegLog <- as.numeric(predict(modelo_RegLog, training_set, type = 'raw'))

AUC_RegLog <- multiclass.roc(training_set$braaksc, prediccion_RegLog)
AUC_RegLog$auc
## Multi-class area under the curve: 0.7797

Para el modelo de Regresión Logística con 8 predictores, devuelve un AUC de 0.7784.



Con este procedimiento obtenemos la siguiente tabla de valores:

Modelo AUC ROC
J48 (8 predictores) 0.8367
J48 (4 predictores) 0.7549
MLP 0.7998
RegLog 0.7784


como el comando multiclass.roc no devuelve réplicas, no podemos hacer tests estadísticos para contrastar si existen diferencias entre nuestros clasificadores, si bien a simple vista parece ser que nuestros modelos tienen un rendimiento similar.


6.4 Contraste de medias para métricas AUC

Por tanto, repetimos el procedimiento, pero esta vez teniendo en cuenta las diversas mediciones realizadas duante el entrenamiento de los modelos, para poder realizar el contraste de medias.


Procedemos al contraste de la normalidad por test de Shapiro-Wilk, cuya hipótesis nula H0 establece que los datos son normales:

apply(metricas_auc, 2, shapiro.test)
## $`J48~AUC`
## 
##  Shapiro-Wilk normality test
## 
## data:  newX[, i]
## W = 0.96873, p-value = 0.8388
## 
## 
## $`J48_2~AUC`
## 
##  Shapiro-Wilk normality test
## 
## data:  newX[, i]
## W = 0.96963, p-value = 0.8525
## 
## 
## $`MLP~AUC`
## 
##  Shapiro-Wilk normality test
## 
## data:  newX[, i]
## W = 0.93938, p-value = 0.3746
## 
## 
## $`RegLog~AUC`
## 
##  Shapiro-Wilk normality test
## 
## data:  newX[, i]
## W = 0.97172, p-value = 0.8826

Si aplicamos la corrección de Bonferroni al alfa, tenemos un valor crítico de 0.05/4 = 0.0125.

Vemos que las métricas AUC de los modelos son normales, ya que todos nuestros p-valores son mayores a 0.0125 y por tanto no podemos rechazar la hipótesis nula.

metricas_auc


A continuación hacemos el contraste para la igualdad de varianzas del AUC. Como son 4 poblaciones normales e independientes, usaremos el test de esfericidad de Bartlett, cuya hipótesis nula H0 establece que se cumple la condición de homocedasticidad en las muestras.

bartlett.test(list(metricas_auc$`J48~AUC`, metricas_auc$`J48_2~AUC`, metricas_auc$`MLP~AUC`, metricas_auc$`RegLog~AUC`))
## 
##  Bartlett test of homogeneity of variances
## 
## data:  list(metricas_auc$`J48~AUC`, metricas_auc$`J48_2~AUC`, metricas_auc$`MLP~AUC`, metricas_auc$`RegLog~AUC`)
## Bartlett's K-squared = 4.7587, df = 3, p-value = 0.1903

Puesto que nuestro p-valor es mayor a 0.05, no rechazamos la hipótesis nula. Vemos por tanto que las métricas AUC de los modelos tienen varianzas iguales.


Puesto que en los datos se cumplen las condiciones de normalidad y homocedasticidad, podemos aplicarles la t de Student o un ANOVA. Por comodidad, haremos uso del ANOVA esta vez.

matriz_anova <- as.data.frame(metricas_auc$`J48~AUC`)
matriz_anova <- cbind(matriz_anova, rep("J48", 15))
colnames(matriz_anova) <- c("AUC","Modelo")


variable_aux <- as.data.frame(metricas_auc$`J48_2~AUC`)
variable_aux <- cbind(variable_aux, rep("J48_2", 15))
colnames(variable_aux) <- c("AUC","Modelo")


variable_aux2 <- as.data.frame(metricas_auc$`MLP~AUC`)
variable_aux2 <- cbind(variable_aux2, rep("MLP", 15))
colnames(variable_aux2) <- c("AUC","Modelo")



variable_aux3 <- as.data.frame(metricas_auc$`RegLog~AUC`)
variable_aux3 <- cbind(variable_aux3, rep("RegLog", 15))
colnames(variable_aux3) <- c("AUC","Modelo")



matriz_anova <- rbind2(matriz_anova, variable_aux)
matriz_anova <- rbind2(matriz_anova, variable_aux2)
matriz_anova <- rbind2(matriz_anova, variable_aux3)
matriz_anova$Modelo <- as.factor(matriz_anova$Modelo)

str(matriz_anova)
## 'data.frame':    60 obs. of  2 variables:
##  $ AUC   : num  0.643 0.709 0.694 0.732 0.77 ...
##  $ Modelo: Factor w/ 4 levels "J48","J48_2",..: 1 1 1 1 1 1 1 1 1 1 ...
anova_auc <- aov(AUC~Modelo, data = matriz_anova)
summary.aov(anova_auc)
##             Df  Sum Sq  Mean Sq F value   Pr(>F)    
## Modelo       3 0.03143 0.010475   9.296 4.39e-05 ***
## Residuals   56 0.06310 0.001127                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


El ANOVA detectó diferencias significativas entre las áreas bajo la curva ROC de los diversos clasificadores. Procedemos con el test de Tukey para ahondar en los resultados del ANOVA:

TukeyHSD(anova_auc)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = AUC ~ Modelo, data = matriz_anova)
## 
## $Modelo
##                    diff          lwr        upr     p adj
## J48_2-J48    0.01706168 -0.015394444 0.04951780 0.5096157
## MLP-J48      0.04251276  0.010056642 0.07496889 0.0054470
## RegLog-J48   0.05951708  0.027060952 0.09197320 0.0000581
## MLP-J48_2    0.02545109 -0.007005037 0.05790721 0.1734197
## RegLog-J48_2 0.04245540  0.009999273 0.07491152 0.0055229
## RegLog-MLP   0.01700431 -0.015451813 0.04946043 0.5124997

Hay diferencias estadísticamente significativas entre el rendimiento de RegLog y J48 (en sus dos variantes) y entre J48 (8 pred.) y MLP (pero curiosamente, la misma comparación entre MLP y J48 con 4 predictores da un p-valor ajustado superior al punto crítico).



7 Últimas palabras

Para finalizar, comentar que quizás podríamos mejorar más aún nuestros modelos si enfocasemos este problema de clasificación como uno de 3 clases, agrupando los estadios braak en los niveles “1-2”, “3-4” y “5-7”. Dicha modificación permitiría mayor flexibilidad a la hora de clasificar el estado de un paciente, ya que los estadios agrupados presentan histopatologías y síntomas cognitivos similares ( Véase Petersen V., Mikkel (2017) ), por no decir que aumentaríamos significativamente la precisión base del problema del 33% actual a un 54,6% (ver cálculo debajo).

C <- length (which(dataset$braaksc=="estadio_5"))
D <- length (which(dataset$braaksc=="estadios_6_7")) 
# Se recuerda que previamente agrupamos a los individuos del estadio 6 y 7 en el nivel 6
precision_actualizada <- ((C+D) / nrow(dataset))*100
print(paste0(precision_actualizada,"%"))
## [1] "54.6456692913386%"