21 min read

Classification and diagnostic prediction of cancers

2018/07/21

Classification and diagnostic prediction of cancers using gene expression profiling and artificial neural networks

En este post se analiza un caso basado en los datos del artículo: (Lantz 2013)

Los datos se pueden obtener directamente de la revista Nature Medicine.

En dicho artículo se investiga la predición del diagnóstico de un tipo de cancer, “small, round blue cell tumors (SRBCTs)” en la infancia usando información del perfil de expresión génica obtenica mediante técnicas de microarrays.

Se estudian 4 tipos de canceres:

The small, round blue cell tumors (SRBCTs) of childhood, which include neuroblastoma (NB), rhabdomyosar- coma (RMS), non- Hodgkin lymphoma (NHL) and the Ewing family of tumors (EWS), are so named because of their similar appearance on routine histology1. However, accurate diagnosis of SRBCTs is essential because the treatment options, responses to therapy and prognoses vary widely depending on the diagnosis. As their name implies, these cancers are difficult to distinguish by light microscopy, and currently no single test can precisely distinguish these cancers.

El proceso de clasificación es mucho más elaborado que los presentados aquí. No reproduciremos este esquema, aunque es importante que se entienda. Se basa en una red neuronal artificial usando 3-fold crossvalidación y repitiendo el proceso 1250 veces mediante particiones aleatorias (proceso similar al bootstrap) para poder estudiar la robustez del modelo. Observar que el número de variables es muy grande (2308) así que se ha optado por realizar un análisis de componentes principales para reducir la dimensión de las variables iniciales y usar solo las 10 primeras en el algoritmo.

El análisis de componentes principales (PCA, en ingles) es una técnica básica y muy utilizada en análisis multivariante para reducir el número de variables creando nuevas variables como combinación lineal de las originales buscando máximizar la varianza explicada. Como no sé si sabeis realizar en R un PCA he optado por crear un fichero con el resultado del PCA llamado pcaComponents.csv.

En este post se usarán los datos del artículo para implementar el algoritmo de red neuronal artificial y “support vector machine” (SVM) para predecir los cuatro tipos de canceres.

Algoritmo Red Neuronal Artificial1

Una Red Neuronal Artificial se basa en establecer la relación entre una serie de señales de entrada (inputs) y de salida (outputs) a partir de un modelo inspirado en el cerebro biológico. Una RNA utiliza una serie de nodos o neuronas para resolver los problemas de aprendizaje. El entrenamiento consiste en las fases de forwardpropagation en el que las neuronas se activan para producir una señal de salida y la fase de backpropagation en el que los valores \(w_i\) que son los pesos de los inputs de cada una de las neuronas se van ajustando para maximizar el alineamiento de la respuesta \(f(X)\) con los valores respuesta del juego de entrenamiento. La no-linealidad se obtiene utilizando las llamadas funciones de activación en cada una de las neuronas. Estas funciones de activación son funciones como la función logística \(a(z)=1/(1+e^{-z})\) o la función ReLU \(a(z)=max(0,z)\).

El ajuste de los pesos se realiza minimizando la función de coste que en general toma la forma de la suma de errores cuadrados \[C(o^{s}_{l})=\frac{1}{2}\sum_{s}\sum_l(t^{s}_{l}-o^{s}_{l})^2\] donde s es el índice de cada uno de los ítems del entrenamiento, o es el valor de salida, t es el valor real de salida y l la neurona de salida. El ajuste se realiza en cada iteración del entrenamiento utilizando el algoritmo de backpropagation: se calcula el gradiente de la función de coste y se ajustan los valores de los \(w_i\), ya que la función de error depende indirectamente de cada uno de los \(w_i\) de la red.

  • PROS

    • Adaptables tanto a problemas de clasificación como de predicción numérica.
    • Capaz de modelar patrones más complejos que cualquier otro algoritmo de ML.
    • Hace pocas suposiciones sobre las relaciones subyacentes de los datos.

  • CONS

    • Computacionalmente costoso y lento de entrenar, especialmente en topologías complejas.
    • Tiende a sobreajustar los datos de entrenamiento.
    • Debido a su no linealidad, se convierte en una “caja negra” que resulta difícil de interpretar.

Algoritmo Support Vector Machine2

Se basa en la generación de unos límites llamados hiperplanos que dividen un espacio de alta/infinita dimensionalidad (dimensiones de los puntos de entrenamiento) en una serie de particiones homogéneas. Estos hiperplanos definen el margen de separación entre las clases (a la vez que definen los límites) e intuitivamente corresponden a la mayor distancia al punto de entrenamiento de cualquier clase.

Una SVM es capaz de ajustar una relación no lineal entre variables ampliando la dimensionalidad utilizando el kernel trick (truco del kernel). De esta manera, no se limita a modelar relaciones lineales entre variables.

  • PROS

    • Adaptables tanto a problemas de clasificación como de predicción numérica.
    • No afectado por datos ruidosos y no tiende a sobreajustar.
    • Conceptualmente más sencillo que las redes neuronales.
    • Ha ganado popularidad debido a su exactitud y a haber ganado grandes competiciones en minería de datos.

  • CONS

    • Encontrar el mejor modelo requiere el testeo con varios kernels y parámetros.
    • Puede lento de entrenar, sobre todo si los datos de entrada tienen un gran número de parámetros o ejemplos.
    • El modelo se convierte en una “caja negra” que resulta difícil de interpretar.

Desarrollo para el clasificador de red neuronal artificial

1. Obtención de datos

Lectura de datos

Leemos los datos a partir de los ficheros en el directorio del documento. El fichero de datos es el fichero data.csv, que no es el resultado del PCA del conjunto de datos, por lo que se deberá realizar antes de la normalización: Casa

# clasificación
tumor.class <- read.csv("class.csv")
# datos de tumores
tumor.data <- read.csv("data.csv")

2. Exploración y preparación de los datos

Principal Component Analysis (PCA)

El PCA consiste en encontrar una serie de componentes independientes entre ellos que sean combinación lineal de las variables. El objetivo es reducir la dimensionalidad de los datos (en este caso, 2308 variables). Los componentes se ordenan en función de la varianza.

Para realizar el PCA, utilizamos la función prcomp del paquete stats en la que podemos indicar el escalado y el centrado (este último parámetro es redundante porque el elcalado ya centra). El escalado significa calcular el valor de las observaciones en unidades de \(\sigma\). A continuación, seleccionamos en la variable pcaData los 10 componentes principales (nuestras nuevas variables) requeridos:

comps <- prcomp(tumor.data, scale. = TRUE, center = TRUE)
pcaData <- comps$x[,1:10]

Podemos hacernos a la idea de cuánta información se va a dejar de utilizar a partir de la proporción acumulada de la varianza para PC10, que es la proporción de varianza explicada respecto a los datos originales. La proporción acumulada es 0.53551 lo que indica que con 10 componentes principales explicamos un 53.551% de la varianza original de los datos. Para hacernos a la idea, con 50 CPs explicaremos alrededor del 90%.

Normalización

Normalizamos los datos para que los datos entre componentes principales sean comparables y las componentes con más varianza no estén sobrerepresentadas:

# Normalización: función
normalize <- function(x) {
  return ((x - min(x)) / (max(x) - min(x)))
}

# pcaData: dataFrame
pcaData <- as.data.frame(pcaData)

# Normalización de los datos
pcaData.norm <- as.data.frame(lapply(pcaData[1:10], normalize))

Codificación

Debemos codificar ortogonalmente la salida de la RN. Establecemos 4 vectores según el tipo de cáncer que queremos clasificar y los asimilaremos a las neuronas de salida en el modelo.Para cada uno de los vectores, uno de las componentes será 1 y el resto 0. Los cánceres tipo EWS, BL, NB, RMS están codificados en las posiciones 1, 2, 3 y 4 respectivamente. Finalmente eliminamos la codificación original de la columna norm.

# tratamiento de etiquetas: codificamos
# EWS, BL, NB, RMS
pcaData.norm <- cbind(pcaData.norm, tumor.class)

pcaData.norm$EWS[pcaData.norm$x == 1] = 1     # EWS
pcaData.norm$EWS[pcaData.norm$x != 1] = 0     # NO EWS
pcaData.norm$BL[pcaData.norm$x == 2] = 1      # BL
pcaData.norm$BL[pcaData.norm$x != 2] = 0      # NO BL
pcaData.norm$NB[pcaData.norm$x == 3] = 1      # NB
pcaData.norm$NB[pcaData.norm$x != 3] = 0      # NO NB
pcaData.norm$RMS[pcaData.norm$x == 4] = 1     # RMS
pcaData.norm$RMS[pcaData.norm$x != 4] = 0     # NO RMS

# borramos
pcaData.norm$x <- NULL

head(pcaData.norm)
##         PC1       PC2       PC3       PC4       PC5       PC6       PC7
## 1 0.5553378 0.6720497 0.5529870 0.7939724 0.4120621 0.5617057 0.3062126
## 2 0.5539219 0.8370178 0.5915564 0.7943331 0.3202740 0.4211253 0.2590429
## 3 0.2953112 0.7405051 0.5635048 0.8039138 0.3314715 0.4644812 0.2024393
## 4 0.4861106 0.9143636 0.6448583 0.7744128 0.4645962 0.8406421 1.0000000
## 5 0.1713988 0.6089361 0.8709542 0.7197806 0.4746517 0.7161501 0.4008877
## 6 0.2326652 0.5098598 0.8711509 0.7618714 0.4947435 0.6124415 0.1275063
##         PC8       PC9      PC10 EWS BL NB RMS
## 1 0.4943006 0.2804608 0.5177177   1  0  0   0
## 2 0.4069738 0.3707326 0.5066594   1  0  0   0
## 3 0.4499325 0.3981021 0.5296861   1  0  0   0
## 4 0.6836750 0.1591436 0.3745376   1  0  0   0
## 5 0.2947759 0.5089563 0.5510486   1  0  0   0
## 6 0.2683417 0.3441790 0.5601478   1  0  0   0

Separación de conjuntos training y test

Para finalizar con el tratamiento de los datos, separamos los conjuntos generando numMuestras índices para el conjunto de train. El conjunto de test será el conjunto complementario de los datos:

# semilla
set.seed(12345)
# buscamos el 67% de los índices (rows)
v <- c(1:nrow(pcaData.norm))
percentage <- 67
numMuestras <- round(nrow(pcaData.norm) * percentage / 100)
# indices de las muestras
index_train <- sample(v, replace = FALSE, size = numMuestras)

Definimos finalmente los juegos de train y test:

train_set <- pcaData.norm[index_train,]
test_set <- pcaData.norm[-index_train,]

3. Entrenamiento del modelo

Se genera el modelo de manera similar, utilizando el parámetro hidden para la red con tres neuronas y con las restricciones del enunciado (una capa). La salida está codificada a partir de los parámetros de tipo:

# una hidden y una neurona
model.1 <- neuralnet(EWS + BL + NB + RMS ~ PC1 + PC2 + PC3 + PC4 + 
                       PC5 + PC6 + PC7 + PC8 + PC9 + PC10, 
                       data = train_set)

plot(model.1, rep = "best", main = "Modelo con una neurona")

# una hidden y tres neuronas
model.3 <- neuralnet(EWS + BL + NB + RMS ~ PC1 + PC2 + PC3 + PC4 + 
                       PC5 + PC6 + PC7 + PC8 + PC9 + PC10, 
                       data = train_set, hidden = 3)

plot(model.3, rep = "best", main = "Modelo con tres neuronas")

4. Evaluación del modelo

Utilizando los datos de test calculamos las predicciones para nuestros modelos.:

# resultado con una neurona en la capa oculta
model_results.1 <- compute(model.1, test_set[1:10])

# resultado con tres neuronas en la capa oculta
model_results.3 <- compute(model.3, test_set[1:10])

A continuación, para facilitar la invocación del método confusionMatrix se recodifican los valores obtenidos. El criterio será el valor de salida con el máximo obtenido. Se han probado otros métodos como el redondeo, pero con este método:

  1. Aseguramos que siempre tengamos una clasificación: utilizando redondeo se daban casos con en el que la salida para un valor de test eran todo 0s.
  2. Evitamos tratar los datos de salida (más allá de lo que lo hacemos) al evitar casos como una neurona de salida a -1.

La función encodeType devuelve la codifiación para el tipo de tumor en función del valor máximo del componente del vector. Los valores de realClass son los valores del conjunto de test real reclasificados (originales) y los valores de clasif1 y clasif3 son los que corresponden a la clasificación de la red neuronal teniendo en cuenta el punto anterior:

encodeType <- function(modelResult) {
  result <- c()
  nrow(result)
  for (v in 1:nrow(modelResult)) {
    y <- modelResult[v,]
    replace(y, y==-1, 0)
    r <- which(max(y) == y)
    result <- append(result, switch (r, "EWS", "BL", "NB", "RMS"))
  }
  return (result)
}

realClass <- as.factor(encodeType(test_set[11:14]))

# resultados de la clasificación en los datos de test
clasif1 <- as.factor(encodeType(model_results.1$net.result))

# resultados de la clasificación en los datos de test
clasif3 <- as.factor(encodeType(model_results.3$net.result))

La matriz de confusión para el clasificador de una neurona:

cm1 <- confusionMatrix(clasif1, realClass)
cm1
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction BL EWS NB RMS
##        BL   0   0  0   0
##        EWS  1   8  0   0
##        NB   2   0  4  12
##        RMS  0   0  0   0
## 
## Overall Statistics
##                                                 
##                Accuracy : 0.4444444             
##                  95% CI : (0.2547988, 0.6467358)
##     No Information Rate : 0.4444444             
##     P-Value [Acc > NIR] : 0.5737806             
##                                                 
##                   Kappa : 0.3076923             
##  Mcnemar's Test P-Value : NA                    
## 
## Statistics by Class:
## 
##                      Class: BL Class: EWS Class: NB Class: RMS
## Sensitivity          0.0000000  1.0000000 1.0000000  0.0000000
## Specificity          1.0000000  0.9473684 0.3913043  1.0000000
## Pos Pred Value             NaN  0.8888889 0.2222222        NaN
## Neg Pred Value       0.8888889  1.0000000 1.0000000  0.5555556
## Prevalence           0.1111111  0.2962963 0.1481481  0.4444444
## Detection Rate       0.0000000  0.2962963 0.1481481  0.0000000
## Detection Prevalence 0.0000000  0.3333333 0.6666667  0.0000000
## Balanced Accuracy    0.5000000  0.9736842 0.6956522  0.5000000

El modelo tiene una modesta precisión (accuracy) de 0.4444444444. Analizando las estadísticas por clase comprobamos que el balance sensibilidad/especificidad es bastante bueno para la clase EWS y en cambio el modelo no detecta ningún paciente de los tipos BL y RMS. Finalmente comentar que el valor de kappa (0.3076923077) es bajo según las tablas de valoración de los índices según la literatura (Lantz 2013).

Se puede concluir que una RNN de una neurona en la capa hidden no es capaz de clasificar correctamente este juego de datos, seguramente porque estamos sobrecargándola en el proceso de entrenamiento: entrenamos primero EWS y no tiene capacidad para ajustar los valores de salida de las otras clases.

Para el clasificador de tres neuronas:

cm3 <- confusionMatrix(clasif3, realClass)
cm3
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction BL EWS NB RMS
##        BL   3   1  0   0
##        EWS  0   6  0   0
##        NB   0   1  4   1
##        RMS  0   0  0  11
## 
## Overall Statistics
##                                                 
##                Accuracy : 0.8888889             
##                  95% CI : (0.7084131, 0.9764725)
##     No Information Rate : 0.4444444             
##     P-Value [Acc > NIR] : 0.000001950483        
##                                                 
##                   Kappa : 0.8421053             
##  Mcnemar's Test P-Value : NA                    
## 
## Statistics by Class:
## 
##                      Class: BL Class: EWS Class: NB Class: RMS
## Sensitivity          1.0000000  0.7500000 1.0000000  0.9166667
## Specificity          0.9583333  1.0000000 0.9130435  1.0000000
## Pos Pred Value       0.7500000  1.0000000 0.6666667  1.0000000
## Neg Pred Value       1.0000000  0.9047619 1.0000000  0.9375000
## Prevalence           0.1111111  0.2962963 0.1481481  0.4444444
## Detection Rate       0.1111111  0.2222222 0.1481481  0.4074074
## Detection Prevalence 0.1481481  0.2222222 0.2222222  0.4074074
## Balanced Accuracy    0.9791667  0.8750000 0.9565217  0.9583333

El modelo tiene una precisión (accuracy) de 0.8888888889. Tras el entrenamiento se observa que han mejorado todos los valores excepto la sensibilidad de EWS. El valor de kappa (0.8421052632) es muy alto.

5. Mejora del rendimiento del modelo

Podemos mejorar el modelo a riesgo de sobreentrenar la red y generar un modelo que se ajuste peor a los datos principalmente:

  • Añadiendo componentes principales a los modelos generados
  • Añadiendo capas/neuronas

Como ejemplo veremos si una red con dos capas ocultas de 3 neuronas cada uno conservando el número de parámetros da mejor resultado:

model.nm <- neuralnet(EWS + BL + NB + RMS ~ PC1 + PC2 + PC3 + PC4 + 
                       PC5 + PC6 + PC7 + PC8 + PC9 + PC10, 
                       data = train_set, hidden = c(3,2) )

model_results.nm <- compute(model.nm, test_set[1:10])
clasif.nm <- as.factor(encodeType(model_results.nm$net.result))
cm.nm <- confusionMatrix(clasif.nm, realClass)

cm.nm
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction BL EWS NB RMS
##        BL   0   0  0   0
##        EWS  1   6  0   1
##        NB   2   2  4   4
##        RMS  0   0  0   7
## 
## Overall Statistics
##                                                 
##                Accuracy : 0.6296296             
##                  95% CI : (0.4236796, 0.8059928)
##     No Information Rate : 0.4444444             
##     P-Value [Acc > NIR] : 0.04101057            
##                                                 
##                   Kappa : 0.4934334             
##  Mcnemar's Test P-Value : NA                    
## 
## Statistics by Class:
## 
##                      Class: BL Class: EWS Class: NB Class: RMS
## Sensitivity          0.0000000  0.7500000 1.0000000  0.5833333
## Specificity          1.0000000  0.8947368 0.6521739  1.0000000
## Pos Pred Value             NaN  0.7500000 0.3333333  1.0000000
## Neg Pred Value       0.8888889  0.8947368 1.0000000  0.7500000
## Prevalence           0.1111111  0.2962963 0.1481481  0.4444444
## Detection Rate       0.0000000  0.2222222 0.1481481  0.2592593
## Detection Prevalence 0.0000000  0.2962963 0.4444444  0.2592593
## Balanced Accuracy    0.5000000  0.8223684 0.8260870  0.7916667

Se observa un empeoramiento del accuracy por lo que se ha sobreentrenado el modelo y los datos de test se ajustan peor que en el caso anterior. De todos modos, la estrategia de mejora pasaría por estudiar la topología de la red y el número de componentes principales.

Modelo de tres nodos con 3-fold crossvalidation

La buena definición de (“Cross-Validation for Predictive Analytics Using R,” n.d.) explica que para evitar que la división de los datos (en este caso 67%/33%) provoque que se generen conjuntos de datos “ruidosos” (sobreajustaríamos el ruido) existen técnicas como la de validación cruzada (cross-validation) en el que el conjunto de datos se particionan en k conjuntos (folds) de aproximadamente el mismo tamaño. Se genera un modelo a partir de todos los folds excepto uno y se calcula un error en la predicción (o un estimador de la bondad de la predicción como el accuracy). Se hace lo mismo para todas las combinaciones de folds y se calcula la media del error (o accuracy). De esta manera podemos estimar el error para diferentes parámetros de modelo.

En nuestro caso, establecemos \(k=3\) pero también consideraremos el parámetro decay \(\lambda\), que es un parámetro de ajuste que penaliza la red neuronal para weight \(w\) altos. El ajuste para cada \(w_i\) durante el entrenamiento es: \[w_i = w_i - \eta(\frac{\partial E}{w_i}-\lambda w_i)\]

Para empezar se recodifican las clases de tumor. A continuación se definen los parámetros del cross-validation utilizando el método trainControl para indicar el método (cv) y el número de folds (3). Después creamos una grid con el número de neuronas (3) y el decay (que toma valores de 0 a 1 con un intervalo de 0.05) que es cada una de las combinaciones de parámetros para las redes que generaremos.

El método train genera los modelos definidos por el frame generado por trainControl y confusionMatrix muestra la media de los samples. Al imprimir el modelo se mostrarán los parámetros óptimos:

## Cross-Validated (3 fold) Confusion Matrix 
## 
## (entries are percentual average cell counts across resamples)
##  
##           Reference
## Prediction   BL  EWS   NB  RMS
##        BL  13.3  0.0  1.2  0.0
##        EWS  0.0 31.3  0.0  0.0
##        NB   0.0  1.2 20.5  1.2
##        RMS  0.0  2.4  0.0 28.9
##                             
##  Accuracy (average) : 0.9398
## Neural Network 
## 
## 83 samples
## 10 predictors
##  4 classes: 'BL', 'EWS', 'NB', 'RMS' 
## 
## No pre-processing
## Resampling: Cross-Validated (3 fold) 
## Summary of sample sizes: 55, 55, 56 
## Resampling results across tuning parameters:
## 
##   decay  Accuracy      Kappa        
##   0.00   0.8311287478  0.76423810189
##   0.05   0.9395943563  0.91665666066
##   0.10   0.9272486772  0.89935047543
##   0.15   0.9149029982  0.88179595653
##   0.20   0.9272486772  0.89857944651
##   0.25   0.8672839506  0.81379950874
##   0.30   0.8077601411  0.72532310643
##   0.35   0.7345679012  0.61320133325
##   0.40   0.7226631393  0.59689133576
##   0.45   0.7226631393  0.59226192367
##   0.50   0.6750440917  0.51883098322
##   0.55   0.6626984127  0.49846018514
##   0.60   0.6146384480  0.42437185999
##   0.65   0.5661375661  0.34636821529
##   0.70   0.4704585538  0.19743114782
##   0.75   0.3866843034  0.06467661692
##   0.80   0.3496472663  0.00000000000
##   0.85   0.3496472663  0.00000000000
##   0.90   0.3496472663  0.00000000000
##   0.95   0.3496472663  0.00000000000
##   1.00   0.3496472663  0.00000000000
## 
## Tuning parameter 'size' was held constant at a value of 3
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were size = 3 and decay = 0.05.

El modelo óptimo corresponde a una red neuronal con decay 0.05 que devuelve un accuracy del 93.9594356261%. Podemos observar el comportamiento de accuracy en función de decay.

plot(model)

Desarrollo para el clasificador de support vector machine

1. Obtención de datos

Borrado de datos en memoria

Para prevenir el uso de frames u otros objetos que hayan podido ser modificados, se borran los valores de las variables en memoria aunque se deban volver a cargar los datos raw, pero manteniendo en memoria los resultados de los modelos cm1cm3ybestModel.nn` que se utilizarán en el resultado final:

# Borramos los valores en memoria.
r <- ls()
r <- r[r!="cm1" & r!="cm3" & r!="bestModel.nn"]
rm(list=r)

Lectura de datos

Al igual que para las RN, leemos los datos a partir de los ficheros en el directorio del documento. El fichero de datos es el fichero data.csv:

# clasificación
tumor.class <- read.csv("class.csv")
# datos de tumores
tumor.data <- read.csv("data.csv")

2. Exploración y preparación de los datos

Separación de conjuntos training y test

# EWS, BL, NB, RMS
tumor.data <- cbind(tumor.data, tumor.class)

tumor.data$x[tumor.data$x == 1] <- "EWS"
tumor.data$x[tumor.data$x == 2] <- "BL"
tumor.data$x[tumor.data$x == 3] <- "NB"
tumor.data$x[tumor.data$x == 4] <- "RMS"

Separamos los conjuntos generando numMuestras índices para el conjunto de train. El conjunto de test será el conjunto complementario de los datos:

# semilla
set.seed(12345)
# buscamos el 67% de los índices (rows)
v <- c(1:nrow(tumor.data))
percentage <- 67
numMuestras.svm <- round(nrow(tumor.data) * percentage / 100)
# indices de las muestras
index_train.svm <- sample(v, replace = FALSE, size = numMuestras.svm)

train_set <- tumor.data[index_train.svm,]
test_set <- tumor.data[-index_train.svm,]

test_vals <- test_set$x
test_vals <- factor(test_vals)
levels(test_vals) <- c("EWS", "BL", "NB", "RMS")

3. Entrenamiento del modelo

Generamos dos modelos: valiladot para el modelo lineal y rbf para radial basis y utilizamos todos los datos del conjunto de entrenamiento sin normalizar y sin reducir la dimensionalidad:

linmodel.svm <- ksvm(x~.,data = train_set, kernel = "vanilladot")
##  Setting default kernel parameters
brfmodel.svm <- ksvm(x~.,data = train_set, kernel = "rbf")

4. Evaluación del modelo

Se debe redondear (función round) para poder comparar los valores. Al hacerlo, obtenemos el hiperplano más próximo para los modelos. Se debe factorizar también para hacer más human friendly el resultado/los niveles de clasificación de los tumores:

linmodel.predict <- predict(linmodel.svm, test_set[,1:2308], type = "response")

# redondeamos y factorizamos la predicción
levels(linmodel.predict) <- c("EWS", "BL", "NB", "RMS")

linmatrix <- confusionMatrix(linmodel.predict, test_vals)
linmatrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction EWS BL NB RMS
##        EWS   3  0  0   0
##        BL    0  8  0   1
##        NB    0  0  4   1
##        RMS   0  0  0  10
## 
## Overall Statistics
##                                                 
##                Accuracy : 0.9259259             
##                  95% CI : (0.7571017, 0.9908999)
##     No Information Rate : 0.4444444             
##     P-Value [Acc > NIR] : 0.0000001806676       
##                                                 
##                   Kappa : 0.8937008             
##  Mcnemar's Test P-Value : NA                    
## 
## Statistics by Class:
## 
##                      Class: EWS Class: BL Class: NB Class: RMS
## Sensitivity           1.0000000 1.0000000 1.0000000  0.8333333
## Specificity           1.0000000 0.9473684 0.9565217  1.0000000
## Pos Pred Value        1.0000000 0.8888889 0.8000000  1.0000000
## Neg Pred Value        1.0000000 1.0000000 1.0000000  0.8823529
## Prevalence            0.1111111 0.2962963 0.1481481  0.4444444
## Detection Rate        0.1111111 0.2962963 0.1481481  0.3703704
## Detection Prevalence  0.1111111 0.3333333 0.1851852  0.3703704
## Balanced Accuracy     1.0000000 0.9736842 0.9782609  0.9166667

Para el modelo RBF:

brfmodel.predict <- predict(brfmodel.svm, test_set[,1:2308], type = "response")

# redondeamos y factorizamos la predicción
levels(brfmodel.predict) <- c("EWS", "BL", "NB", "RMS")

brfmatrix <- confusionMatrix(brfmodel.predict, test_vals)
brfmatrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction EWS BL NB RMS
##        EWS   3  0  0   0
##        BL    0  8  0   5
##        NB    0  0  4   3
##        RMS   0  0  0   4
## 
## Overall Statistics
##                                                 
##                Accuracy : 0.7037037             
##                  95% CI : (0.4981863, 0.8624734)
##     No Information Rate : 0.4444444             
##     P-Value [Acc > NIR] : 0.005852564           
##                                                 
##                   Kappa : 0.6                   
##  Mcnemar's Test P-Value : NA                    
## 
## Statistics by Class:
## 
##                      Class: EWS Class: BL Class: NB Class: RMS
## Sensitivity           1.0000000 1.0000000 1.0000000  0.3333333
## Specificity           1.0000000 0.7368421 0.8695652  1.0000000
## Pos Pred Value        1.0000000 0.6153846 0.5714286  1.0000000
## Neg Pred Value        1.0000000 1.0000000 1.0000000  0.6521739
## Prevalence            0.1111111 0.2962963 0.1481481  0.4444444
## Detection Rate        0.1111111 0.2962963 0.1481481  0.1481481
## Detection Prevalence  0.1111111 0.4814815 0.2592593  0.1481481
## Balanced Accuracy     1.0000000 0.8684211 0.9347826  0.6666667

Se observa que el ajuste lineal de la SVM con un accuracy de 0.4444444444 mejora el resultado de la función BRF obtiene peores resultados que el ajuste lineal (0.8888888889). En ambos modelos (lineal/rbf) el valor de la sensibilidad es de 1 en todas las clases excepto en RMS. De todos modos, el valor de la especificidad indica que el único en el que \(sensibilidad=especificidad=1\) es EWS.

5. Mejora del rendimiento del modelo

Para mejorar el modelo podemos estudiar la posibilidad de variar el kernel de la SVM, aunque en este caso el ajuste lineal da muy buenos resultados.

Modelo de tres nodos con 3-fold crossvalidation

Los comentarios sobre el funcionamiento del método de cross-validation del apartado para la red neuronal para obtener los parámetros y estimar el accuracy de los modelos son válidos también para el caso que estemos estudiando el funcionamiento de una SVM. En este se parametrizará la SVM (más allá del kernel único) por lo que solamente habrá un resultado.

El valor del tipo de tumor ya está factorizado del apartado del crossvalidation de la RNN: se usan los mismos datos de entrada.

## Cross-Validated (3 fold) Confusion Matrix 
## 
## (entries are percentual average cell counts across resamples)
##  
##           Reference
## Prediction   BL  EWS   NB  RMS
##        BL  13.3  0.0  0.0  0.0
##        EWS  0.0 32.5  0.0  0.0
##        NB   0.0  0.0 21.7  1.2
##        RMS  0.0  2.4  0.0 28.9
##                             
##  Accuracy (average) : 0.9639
## Support Vector Machines with Linear Kernel 
## 
##   83 samples
## 2308 predictors
##    4 classes: 'BL', 'EWS', 'NB', 'RMS' 
## 
## No pre-processing
## Resampling: Cross-Validated (3 fold) 
## Summary of sample sizes: 55, 55, 56 
## Resampling results:
## 
##   Accuracy      Kappa       
##   0.9638447972  0.9500100827
## 
## Tuning parameter 'C' was held constant at a value of 1

El modelo óptimo corresponde a una SVM que devuelve un accuracy del 96.3844797178%.

Resultados finales y elección de modelo

Podemos resumir los resultados de los modelos obtenidos en la siguiente tabla:

Método Parámetros Accuracy
Red Neuronal hidden 1 neurona 0.4444444444
Red Neuronal hidden 3 neuronas 0.8888888889
Red Neuronal cross-validation 0.9395943563
SVM Lineal 0.9259259259
SVM BRF 0.7037037037
SVM cross-validation (Lineal) 0.9638447972

La red neuronal con 1 neurona en capa oculta obtiene resultados muy pobres tal como se ha comentado, pero la validación cruzada indica que existen configuraciones que optimizan el resultado de la red con 3 neuronas en capa oculta (variando el parámetro decay). El modelo de SVM con kernel lineal es el mejor entre los modelos no obtenidos a partir de cross-validation y el cv de la SVM lineal devuelve el mejor accuracy global. Se concluye así que los modelos con una mayor accuracy son las SVM lineales, por lo que sería el método escogido.

Referencias

“Cross-Validation for Predictive Analytics Using R.” n.d. http://www.milanor.net/blog/cross-validation-for-predictive-analytics-using-r/.

Lantz, Brett. 2013. Machine Learning with R. Packt Publishing Ltd.

  1. A partir de la definición de (Lantz 2013) y de otras fuentes online (Wikipedia, blogs, etc.)↩

  2. A partir de la definición de (Lantz 2013) y de otras fuentes online (Wikipedia, blogs, etc.)↩