---
title: "Clasificación IV - Software R"
author:
  - Francisco Charte Ojeda - <a href="https://fcharte.com">Sitio web</a>
  - Universidad de Jaén - DaSCI
date: "`r format(Sys.time(), '%d %B %Y')`"
output:
  html_notebook:
    toc: yes
subtitle: Seminario Permanente de Formación en Inteligencia Artificial Aplicada a
  la Defensa
---

# Introducción 

Este es un _notebook_ R creado con RStudio, mediante la opción **File>New File>R Notebook**. Se trata de un tipo de documento que permite combinar texto y formato con código en lenguaje R y los resultados que produce dicho código.

Aunque en el presente _notebook_ se recurrirá a representaciones a gráficas para llevar a cabo un análisis exploratorio de los datos cuando sea preciso, los detalles sobre cómo confeccionar gráficas y obtener indicadores estadísticos de los datos **se facilitan en un documento separado**. En este se describen únicamente los mecanismos para construir, evaluar y usar diferentes tipos de clasificadores.

## Paquetes R para clasificación

En CRAN existe un _Task View_ específico sobre [aprendizaje automático](https://cran.r-project.org/web/views/MachineLearning.html) con un centenar de paquetes que implementan una gran parte de los algoritmos de clasificación propuestos en la literatura, desde los más clásicos como los árboles de decisión o kNN a métodos basados _ensembles_ y aprendizaje profundo.

Por regla general, todos los paquetes ofrecen un método `train()`, encargado de ejecutar el algoritmo de aprendizaje sobre las muestras facilitadas y devolver el correspondiente modelo, y un método `test()` que, a partir de dicho modelo, predice la clase para nuevas instancias.

# Clasificación de conjuntos de datos muy simples

Supongamos que nuestros datos tienen una estructura muy sencilla, contando con un par de variables predictoras, además del atributo que actúa como etiqueta de clase, y que se trata de un problema de clasificación binario. Podemos construir este conjunto de datos a partir de `iris`, quedándonos con dos variables y solo dos de las tres clases de flor:

```{r}
data(iris)
set.seed(73)

# Tomamos tres de los cinco atributos y las muestras que no sean de la familia `setosa`
datos <- iris[iris$Species != "setosa", c("Petal.Length", "Petal.Width", "Species")]

# Convertimos el atributo de clase en un número
datos$Species <- as.numeric(datos$Species) - 1

datos[sample(nrow(datos), 10), ] # Vistazo a 10 muestras tomadas aleatoriamente
```

Aunque tenemos muy pocos datos, solo 100 muestras para entrenar y evaluar el clasificador, vamos a usar el enfoque de particionamiento _hold-out_ por simplicidad:

```{r}
set.seed(7)
test.ind <- sample(nrow(datos), size = nrow(datos) * 0.25)
train <- datos[-test.ind, ]
test <- datos[test.ind, ]

str(train)
str(test)
```

Por último, vamos a examinar la distribución de las muestras de test recurriendo a una gráfica de dispersión:

```{r}
library(ggplot2)

ggplot(train, aes(x = Petal.Length, y = Petal.Width)) +
  geom_point(aes(color = as.factor(Species))) +
  stat_smooth(method = lm, color= "black", se = FALSE)
```

Como puede apreciarse, las muestras pertenecientes a cada familia podrían separarse casi perfectamente si usáramos como frontera una línea perpendicular a la mostrada y que la cruzase aproximadamennte por el punto central. Esto denota que un clasificador lineal simple, creado a partir de la línea regresión, sería suficiente.

## Clasificación mediante regresión lineal

Vamos a ajustar un regresor lineal con el objeto de predecir la especie, que es un valor numérico en nuestro conjunto de datos, usando el resto de variables disponibles:

```{r}
modelo <- lm(Species ~ ., data = train)
modelo
```

El modelo resultando nos facilita los coeficientes de regresión para cada variable, así como el _intercept_ o punto de corte de la recta con el eje Y o de ordenadas. Estos coeficientes son accesibles individualmente como `modelo$coefficients[N]`. Esto nos permite componer la fórmula de predicción y, con ella, construir nuestro clasificador:

```{r}
prediccion <- 
  modelo$coefficients[1] + 
  modelo$coefficients[2] * test$Petal.Length + 
  modelo$coefficients[3] * test$Petal.Width

test$Species
prediccion
round(prediccion)
```

Como puede apreciarse, la predicción es un valor real dado que estamos usando una fórmula de regresión. No obstante, podemos convertirla en una etiqueta de clase simplemente redondeando al entero más próximo. Esto sería lo mismo que decir que aplicamos un umbral de corte, de forma que los valores inferiores a `1.5` se consideran de la clase 1 y los mayores a ese valor de clase 2.

Comparando los valores predichos con los reales podemos calcular el ratio de acierto _Accuracy_ que, como se aprecia a continuación, es perfecto sobre los datos de test. Esto se debe a la sencillez del conjunto de datos original.

```{r}
accuracy <- sum(test$Species == round(prediccion)) / nrow(test) * 100 
paste0(accuracy, "% de acierto")
```

## Clasificación mediante agrupamiento

Las técnicas de _clustering_ son no supervisadas, es decir, no precisan de una etiqueta de clase para asignar cada muestra del conjunto de datos a un grupo u otro, guiándose exclusivamente por las distancias entre los puntos en el espacio n-dimensional generado por las variables que describen cada instancia.

Por ello, usamos solo las dos variables que indican la longitud y anchura de pétal para ejecutar K-medias, el algoritnmo de agrupamiento más conocido y simple:

```{r}
library(cluster)

set.seed(73)
clusters <- kmeans(train[,1:2], 2)

clusters
```

En la información devuelta observamos las coordenadas que actúan como centro de los dos _cluster_ generados, la asignación de las muestras a uno u otro y la suma de error cuadrático interno que indica el grado de similitud de los puntos. Podemos representar los datos para apreciar el trabajo hecho por K-medias.

```{r}
ggplot(train, aes(x = Petal.Length, y = Petal.Width)) +
  geom_point(aes(color = as.factor(clusters$cluster), shape = as.factor(Species)))

clusplot(train[,1:2], clusters$cluster, color = TRUE)
```

En la nube de puntos se aprecia cómo en la parte central, donde los dos _cluster_ se entremezclan, tenemos puntos de diferente color (que denota el grupo al que pertenece) y forma (que denota la especie, no usada para hacer el agrupamiento). La segunda representación gráfica nos permite apreciar las fronteras que separan a un grupo de otro.

Si bien la finalidad del agrupamiento es más descriptivo que predictivo, podemos usar los centroides de los _cluster_ como referencia, calculando la distancia de las nuevas muestras (las de test) a cad centro y quedándonos con la menor como etiqueta de clase. Para ello comenzamos escribiendo una función que, recibiendo comoparémtros las muestras de test y los _cluster_, devuelve como resultado la predicción:

```{r}
predict.kmeans <- function(muestras, clusters) {
  diferencias <- t(apply(muestras, 1, function(f) 
    sqrt(rowSums((f - clusters$centers)^2))))
  ifelse(diferencias[,1] < diferencias[,2], 1, 2)
}

```

Usando la anterior función obtenemos las predicciones tanto para las muestras de entrenamiento como las de test y mostrando el ratio de acierto obtenido.

```{r}
train$Prediccion <- predict.kmeans(train[,1:2], clusters)
test$Prediccion <- predict.kmeans(test[,1:2], clusters)

accuracy <- sum(train$Species == train$Prediccion) / nrow(train) * 100 
paste0(accuracy, "% de acierto - entrenamiento")
accuracy <- sum(test$Species == test$Prediccion) / nrow(test) * 100 
paste0(accuracy, "% de acierto - test")
```

# Predicción con árboles de decisión

Los árboles de decisión son modelos de clasificación altamente interpretables, a partir de los cuales es fácil extraer un conjunto de reglas que una persona puede aplicar directamente. En R hay disponibles múltiples algoritos para la generación de árboles, entre ellos el paquete `rpart`. El objetivo será predecir qué pasajeros del Titanic sobrevivirán a partir de datos como su sexo, edad o la clase del pasaje con el que viajaban.

## El conjunto de datos `titanic`

Comencemos explorando el conjunto de datos con el que vamos a trabajar, alojado en la variable `titanic_train` del paquete `titanic`. El primer paso será examinar las variables que tiene dicho conjunto de datos:

```{r}
library(titanic)
datos <- titanic_train

str(datos)
```

De los atributos con que cuenta cada muestra hay varios que no tienen valor como variables predictoras. Es el caso de `PassengerId`, que identifica de forma única a cada muestra; `Name`, con el nombre de los pasajeros; `Ticket`, `Cabin` o `Embarked`. Las variables que nos serán útiles son las siguientes:

- `Pclass`: 1, 2 o 3 indicando la clase del pasaje 
- `Sex`: valores `female` y `male` indicando el sexo
- `Age`: edad de cada pasajero
- `Fare`: precio del pasaje
- `SibSp`: número de hermanos, esposo/esposa a bordo
- `Parch`: número de padres/hijos a bordo
- `Survived`: 0 si el pasajero no sobrevivió, 1 en caso contrario

Nos quedamos solo con esas variables:

```{r}
variables <- c("Pclass", "Sex", "Age", "Fare", "SibSp", "Parch", "Survived")

datos <- datos[, variables]

str(datos)
```

## Análisis de correlaciones entre variables

Dado que tenemos siete variables, repersentar estos datos como puntos en un plano no es fácil. No obstante, podemos analizar las correlaciones existentes entre ellas de una forma mucho más cómoda, mediante un diagrama de correlaciones:

```{r}
library(dplyr)
library(corrplot)

datos %>%
  na.omit %>%
  mutate(Sex = as.numeric(as.factor(Sex))) %>%
  cor %>%
  corrplot(method = "square", type = "upper")
```

Fijándonos en la última columna, que es la actúa como clase de las muestras, observamos que hay una alta correlación negativa con la variable `Sex` y algo menor con la variable `Pclass`. También existe un grado de correlación importante, en este caso positiva, respecto a la variable `Fare`.

## Generación del árbol y evaluación 

Cargamos el paquete `rpart`, dividimos los datos en un subconjunto de entrenamiento y otro de test, generamos el árbol de decisión y mostramos información sobre su estructura:

```{r}
library(rpart)

set.seed(73)
test.ind <- sample(nrow(datos), size = nrow(datos) * 0.25)
train <- datos[-test.ind, ]
test <- datos[test.ind, ]

modelo <- rpart(Survived ~ ., data = train, method = "class")
summary(modelo)
```

En este resumen se aprecian aspectos como la importancia relativa de las variables a la hora de tomar decisiones, el número de muestras en cada nodo, reparto por clases y otros parámetros. En general, resulta mucho más sencillo analizar esa información si la representamos gráficamente tal y como se muestra a continuación:

```{r}
library(rpart.plot)
rpart.plot(modelo)
```

A continuación evaluamos el clasificador usando las muestras de test que habíamos reservado anteriormente. La matriz de confusión nos permite calcular fácilmente diferentes métricas de evaluación, entre ellas _Accuracy_.

```{r}
prediccion <- predict(modelo, test, type = "class")
matriz.confusion <- table(test$Survived, prediccion)
accuracy <- sum(diag(matriz.confusion)) / sum(matriz.confusion)

matriz.confusion
accuracy
```

## Extracción de las reglas de clasificación

Aunque el modelo devuelto por `rpart()` puede emplearse directamente para obtener predicciones futuras, simplemente facilitando los valores de las variables predictoras, puede interesarnos tener un conjunto de reglas que haga posible aplicar el modelo manualmente. Con esa idea, usamos el paquete `tidyrules`  a fin de extraer las reglas y generar una tabla con ellas:

```{r}
library(tidyrules)
library(pander)

reglas <- tidyRules(modelo)
pandoc.table(reglas)
```
## Parámetros de control del árbol

En general, al ejecutar un algoritmo que construye un modelo a partir de datos de ejemplo podemos controlar ciertos aspectos alterando los valores por defecto para ciertos parámetros. En el caso de los árboles de decisión, por ejemplo, podemos modificar el número mínimo de muestras que deberá contener cada nodo hoja, la profundidad máxima del árbol, etc. En el caso de `rpart()` esos parámetros se facilitan como una lista `rpart.control`, tal y como se muestra en el siguiente ejemplo:

```{r}
modelo <- rpart(Survived ~., data = train, method = "class",
                 control = rpart.control(
                   cp = 0.001, 
                   minbucket = 10, 
                   maxdepth = 3))
rpart.plot(modelo)

prediccion <- predict(modelo, test, type = "class")
matriz.confusion <- table(test$Survived, prediccion)
accuracy <- sum(diag(matriz.confusion)) / sum(matriz.confusion)

matriz.confusion
accuracy

reglas <- tidyRules(modelo)
pandoc.table(reglas)
```

Como puede comprobarse, el árbol es mucho más sencillo que el obtenido anteriormente. Su ratio de acierto ha descendido ligeramente, pero a cambio se tiene un modelo más fácilmente aplicable y que también genera un conjunto de reglas más reducido, lo cual también facilita su aplicación.

# Predicción con kNN y estimación realista del rendimiento

La predicción de etiquetas para nuevas muestras de datos basándose en las que corresponden a sus vecinos más cercanos, aquellas instancias ya conocidas y cuya _cercanía_ podemos determinar usando alguna medida de distancia, representa un enfoque fácil de entender. El algoritmo para llevar a cabo esta tarea, conocido como kNN, no construye en general un modelo a partir de los datos de entrenamiento, por lo que todo el trabajo recae en la fase de clasificación. En esta sección usaremos este tipo de clasificador para determinar el tipo de cristal a partir de un conjunto de características químicas. Además, también se explicará cómo llevar a cabo una estimación realista del rendimiento de este clasificador.

## Carga y análisis exploratorio

El conjunto de datos `glass` está disponible en el repositorio UCI, por lo que podemos descargarlo directamente (debemos contar con una conexión a Internet) y alojarlo en memoria. Establecemos el nombre de las variables y procedemos con la fase de análisis exploratorio:

```{r}
glass <- read.csv("http://archive.ics.uci.edu/ml/machine-learning-databases/glass/glass.data")
colnames(glass) <- c("Id", "RI","Na","Mg","Al","Si","K","Ca","Ba","Fe","Class")

str(glass)
```

La variable `Id` es un identificador único para cada muestra de datos, por lo que no resulta de interés para realizar predicciones. La variable `Class` es la que determina la clase, el tipo de cristal, y contempla siete etiquetas distintas. Estamos, por tanto, ante un problema multiclase.

Quedándonos solo con las variables que nos interesan, únicamente vamos a usar un histograma por variable para completar el análisis exploratorio:

```{r}
library (tidyverse)

variables <- c( "RI","Na","Mg","Al","Si","K","Ca","Ba","Fe","Class")
datos <- glass[, variables]

datos %>% 
  gather() %>% 
  ggplot(aes(value)) +
    facet_wrap(~ key, scales = "free") +
    geom_histogram()
```

Lo primero que podemos apreciar es que prácticamente ninguna de las variables sigue una distribución normal en sus valores. Por otra parte, la variable `Class` denota que hay un cierto desequilibro entre clases: se trata de un problema de **clasificación desbalanceada**.

## Parámetros para controlar el particionamiento y evaluación

Aunque con R podríamos llevar a cabo de forma manual el particionamiento de los datos empleando las diferentes técnicas ya conocidas, como _hold-out_, validación cruzada o _leave-one-out_, esta es una tarea ya automatizada en el paquete `caret`. En él, además, también encontramos implementaciones de la mayoría de algoritmos de aprendizaje más populares, incluyendo kNN, por lo que podempos directamente entrenar y evaluar el predictor para obtener una estimación de su rendimiento:

```{r}
library(caret)

glass$Class <- as.factor(glass$Class)

train_control <- trainControl(method = "cv", number = 5, savePredictions = TRUE)

modelo.knn <- train(Class ~ ., data = glass, method = 'kknn',
                  trControl = train_control)
                  
modelo.knn$results

summary(modelo.knn)
```

Usando 5fcv (_5 folds cross validation_) se divide el conjunto de datos en cinco partes de igual número de muestras. En cada iteración se emplearán cuatro de esas partes como entrenamiento y una como test, obteniendo cinco evaluaciones distintas del clasificador. Los datos de rendimiento obtenidos son valores promedio de esas cinco evaluaciones, razón por la que también se incluye la desviación estándar.

En el resumen de resultados se informa de los valores que se han dado a algunos parámetros, como es el caso de `kmax` que determina el número de vecinos tomados en cuenta para realizar la predicción. En la salida resumida se concluye que `k = 5` es la mejor opción. No se han probado otras alternativas para este parámetro y otros potencialmente interesantes, como el tipo de distancia o la función _kernel_ a usar para ponderar las distancias de los vecinos. 

## Ajuste de parámetros del algoritmo

Además de llevar a cabo la estimación del rendimiento mediante la técnica de validación que nos interese, con `caret` también podemos llevar a cabo un proceso de validación interno para la selección de los mejores valores a emplear pra el clasificador. En el caso de kNN, por ejemplo, podemos probar distintos números de vecinos, funciones de distancia y _kernels_  para los pesos:

```{r}
parametros.knn <- expand.grid(kmax = c(1:5), distance = c(1, 2), 
                        kernel = c('gaussian', 'epanechnikov', 'rectangular'))

modelo.knn <- train(Class ~ ., data = glass, method = 'kknn',
                  trControl = train_control, tuneGrid = parametros.knn)
                  
modelo.knn$results

summary(modelo.knn)

```

El ajuste de estos parámetros nos permite elevar el nivel de acierto del clasificador desde el algo más del 86% previo a más del 93%.

# Predicción con máquinas de vectores soporte

Las conocidas como SVM son una familia de modelos predictivos que, si bien resultan menos interpretables que las opciones previas, suelen destacar por su gran rendimiento. Originalmente se diseñaron para abordar problemas de clasificación binaria, por lo que a continuación veremos cómo aplicarlas con el objetivo de identificar automáticamente el correo basura o _spam_.

Comenzamos, como es habitual, llevando a cabo un proceso de análisis exploratorio de los datos para, a continuación, crear un modelo de clasificación que en este caso será una SVM.

## Carga de datos y análisis exploratorio

Desde el [repositorio UCI](http://archive.ics.uci.edu/ml/datasets/spambase) podemos descargar un archivo conteniendo un _dataset_ correspondiente a varios miles de correos electrónicos clasificados como _spam_ y no _spam_. Se trata de un conjunto con 4600 muestras recopilado en Hewlett-Packard Labs. En R podemos importar este mismo _dataset_ desde el paquete `kernlab`:

```{r}
library(kernlab)
data(spam)
```

En la documentación en línea encontraremos información sobre la procedencia del conjunto de datos. También podemos echar un vistazo a sus primeras muestras, así como obtener un resumen de la distribución estadística de los valores que toman sus variables. De esta forma componemos una visión general de la estructura del _dataset_:

```{r}
?spam   # Documentación integrada sobre el dataset

head(spam)
summary(spam)
```

A continuación recurrimos a distintas representaciones gráficas con el objetivo de conocer mejor ciertos aspectos, como la distribución de las clases, frecuencia de ciertos símbolos y palabras en los correos electrónicos según su tipo:

### Cantidad de correos según tipo: spam/legítimo

Comenzamos comprobando la cantidad de correos que hay en el conjunto de datos de cada tipo, para lo cual es suficiente con una gráfica de barras:

```{r}
library(ggplot2)

ggplot(spam, aes(x = type)) + geom_bar(aes(fill = type))
```

### Calculamos la ratio de desbalanceo

Es fácil apreciar que existe un cierto desequilibrio en el número de muestras de cada clase. Calculamos el conocido com IR (_Imbalance Ratio_) para determinar exactamente la medida en que el _dataset_ está desbalanceado:

```{r}
library(dplyr)
conteo <- spam %>% group_by(type) %>% count
conteo$n[1] / conteo$n[2]
```

### Frecuencias de mayúsculas y uso del símbolo $

En los correos basura es habitual abusar de las mayúsculas, así como del uso de ciertos símbolos. Podemos comprobarlo examinando la frecuencia de algunos de los atributos del conjunto de datos diferenciando por tipo de correo:

```{r}
ggplot(spam, aes(x = capitalLong)) + 
  geom_histogram(aes(fill = type), binwidth = 20) +
  scale_x_continuous(limits = c(0, 1750)) +
  scale_y_continuous(limits = c(0, 350))
```

Aquí comprobamos cómo las secuencias largas de texto en mayúsculas corresponden casi por completo a correos basura, mientas que los legítimos destacan por un uso moderado. 

Veamos qué ocurre con el uso del símbolo _$_:

```{r}
ggplot(spam, aes(y = charDollar)) + 
  geom_boxplot(aes(color = type))

spam %>% group_by(type) %>% summarize(Apariciones = sum(charDollar), Promedio = mean(charDollar))
```

En los diagramas de cajas y bigotes se aprecia que su uso es mucho más frecuente en los correos basura, algo qe confirman los valores absolutos y relativos de apariciones calculados a continuación.

### Frecuencias de algunas palabras según el tipo de correo

Es habitual que en los correos basura aparezcan ciertas palabras que son mucho menos corrientes en mensajes legítimos. Vamos a tomar parte de los términos del dataset, concretamente los primeros veinte, y a representarlos en gráficas de barras según el tipo de correo:

```{r}
library(reshape2)

datos <- melt(spam, measure.vars = 1:20, id.vars = "type")
levels(datos$type) <- c("L", "B")
str(datos)

ggplot(datos, aes(y = value, x = type, fill = type)) +
  geom_bar(stat = "identity") +
  facet_wrap(~variable, scales = "free")
```

Observando estas gráficas puede apreciarse que términos como _credit_, _business_ o _receive_ apenas aparecen en los mensajes legítimos pero sí en el correo basura. Lo inverso ocurre con otras palabras (aquí solo se han examinado algunas de las que aparecen en cada mensaje). 

## Introducción al paquete `caret`

Una vez que conocemos el conjunto de datos con el que estamos trabajando, procedemos a entrenar una SVM con la finalidad de generar un modelo de clasificación adecuado. 

En este caso vamos a usar el paquete `caret`. Al igual que el `scikit-learn` de Python, este paquete R nos ofrece una metodología común para usar prácticamente cualquier algoritmo de aprendizaje automático. Para ello se apoya fundamentalmente en dos funciones:

- `train()`: se encarga de entrenar el modelo indicado por el parámetro `method` con el objetivo de crear un predictor acorde a la fórmula entregada como primer argumento. Es capaz de llevar a cabo el particionamiento y estimación del rendimiento del predictor, todo ello de forma automática. Actualmente `caret` nos ofrece más de 230 métodos, cuyos nombres es fácil obtener:

```{r}
library(caret)
names(getModelInfo())
```

Además de particionar los datos, entrenar y evaluar los modelos, el método `train()` también es capaz de ajustar los parámetros de funcionamiento de los algoritmos siempre que se le facilite, mediante el parámetro `tuneGrid`, la lista de parámetros e intervalos a probar. Asimismo, `train()` puede aplicar distintos preprocesamiento a los datos antes de usarlo, por ejemplo normalizando los valores de las variables.

- `predict()`: una vez que tenemos el modelo entrenado, la función `predict()` lo usará para, aplicándolo sobre nuevas muestras de datos, generar las predicciones correspondientes. 

## Clasificación con SVM

Procedemos particionando nuestro conjunto de datos. Para ello no realizaremos un particionamiento aleatorio, como en ejemplos previos, sino de tipo estratificado, de forma que la distribución de clases sea similar entre los datos de entrenamiento y prueba. De esto se encarga la función `createDataPartition()`. Nos reservamos un 30% de los correos para verificar al final el rendimiento del clasificador:

```{r}
library(caret)
set.seed(73)

# Particionamos los datos de forma estratificada, reservando un 30% para test final
indices.test <- createDataPartition(y = spam$type, p = 0.3, list = FALSE)
entrenamiento <- spam[-indices.test, ]
test <- spam[indices.test, ]

conteo.entrenamiento <- entrenamiento %>% group_by(type) %>% count
conteo.test <- test %>% group_by(type) %>% count
cat("Proporción de clases en entrenamiento:", conteo.entrenamiento$n[1] / conteo.entrenamiento$n[2], "\n")
cat("Proporción de clases en test:", conteo.test$n[1] / conteo.test$n[2])
```

A continuación, usando solo los datos de entrenamiento, llevamos a cabo una validación cruzada con 5 porciones a fin de entrenar la SVM y obtener una estimación inicial de su rendimiento. Usamos un _kernel_ de tipo RBF para la proyección en dimensiones superiores y optimizamos la métrica ROC durante el entrenamiento:

```{r}
# Entrenamos un modelo de SVM con kernel RBF y obtenemos el estimador
train_control <- trainControl(
  method = "cv", number = 5, 
  classProbs = TRUE, savePredictions = TRUE,
  summaryFunction = twoClassSummary)

svm.spam <- train(type ~., data = entrenamiento, method = "svmRadial",
      trControl = train_control,
      metric = "ROC")
svm.spam
```

La salida del clasificador nos indica que la mejor configuración es con `C = 1`, alcanzando un ROC de casi el 97%. Podemos obtener métricas adicionales de rendimeinto con el método `multiClassSummary()`:

```{r}
multiClassSummary(svm.spam$pred, svm.spam$levels)
```


Podemos examinar cómo ha cambiado el rendimiento del clasificador en proporción al coste (parámetro C):
```{r}
plot(svm.spam)

```

También annalizar la curva ROC del clasificador generado:
```{r}
library(MLeval)
evalua <- evalm(svm.spam, showplots = FALSE)
evalua$roc
```


Asumiendo que tuviésemos nuevos correos entrantes, usaríamos la función `predict` para determinar si son _spam_ o no. La matriz de confusión nos indica el nivel de acierto con que se realiza esa predicción usando muestras de datos antes no vistas:

```{r}
predicciones <- predict(svm.spam, test, type = "prob")
predicciones$obs <- test$type
predicciones$pred <- predict(svm.spam, test)
multiClassSummary(predicciones, lev = levels(predicciones$obs))

confusionMatrix(predicciones$pred, test$type)
```

Mediante una gráfica de nube de puntos podemos mostrar gráficamente el volumen de aciertos y fallos para cada clase. Básicamente es una visualización de la mamtriz de confusión:

```{r}
pred <- data.frame(pred = predicciones$pred, obs = predicciones$obs)
ggplot(pred, aes(x = pred, y = obs)) + 
  geom_jitter(position = position_jitter(width = 0.25, height = 0.25))
```


# Ensembles de clasificadores

En sesiones previas conocimos distintos tipos de _ensembles_, conjuntos de múltiples clasificadores, entrenados con subconjuntos de muestras y/o variables para inducir diversidad, que combinan sus salidas mediante alguna estrategia de voto. 

En esta última sección, usando el paquete `caret` que ya conocemos, se introducen varios de los ensembles disponibles en R.

## Preparación de paquetes y datos

Comenzamos cargando los paquetes necesarios y particionamos los datos reservando un 30% de muestras para el test final, simulando que tenemos nuevos correos a clasificar:

```{r}
library(caret)
library(randomForest)
trellis.par.set(caretTheme())  # Aspecto para algunas gráficas

# Particionamos los datos
set.seed(73)
indices <- createDataPartition(spam$type, p = .70, list = FALSE)
training = spam[indices, ]
test = spam[-indices, ] 
```

## Entrenamos los modelos multi-clasificadores

Usando en este caso una 10fcv, entrenamos varios _ensembles_ usando los enfoques de _bagging_, _boosting_  y _random forest_ que conocimos en una sesión previa, **Aviso**: este proceso puede tardar un cierto tiempo, al contar con varias miles de muestras que han de usarse para entrenar cientos de modelos distintos según el tipo de ensemble:

```{r}
# Configuramos la selección del modelo durante entrenamiento
train10CV <- trainControl(method = "cv", number = 10, classProbs = TRUE)

# Bagging
bag <- train(type ~ ., data = training, method = "treebag", trControl = train10CV)     
# Boosting
gbm <- train(type ~ ., data = training, method = "gbm", trControl = train10CV, verbose = FALSE) 
# Random forest
rf <- train(type ~ ., data = training, method = "rf", trControl = train10CV)           
```

## Examinamos algunos aspectos del entrenamiento

Obtenemos un resumen de los modelos entrenados y una visualización de algunos detalles. Comenzamos con el _bagging_:

```{r}
bag
```

La estimación nos dice que alcanza casi un 93% de acierto. Veamos en el caso de _random forest_:

```{r}
rf
```

En este caso la ratio de acierto sube a algo más del 94%. Este nivel se alcanza cuando `mtry = 2`, descendiendo para valores superiores. Dicho parámetro determina el número de variables que de forma aleatoria se prueba como candidatas para ir dividiendo los árboles.

Este modelo nos permite, además, observar detalles como el cambio en rendimiento según el número de predictores seleccionados (el citado parámetro `mtry`) o el número de árboles empleados en el _ensemble_:

```{r}
plot(rf)

leyenda <- attr(rf$finalModel$err.rate, "dimnames")[[2]]
plot(rf$finalModel, main = "Error vs Número de árboles")
legend("top", legend = leyenda, lty = 1:3, col = 1:3)
```

Examinemos ahora la estimación de rendimiento para el _boosting_:

```{r}
gbm
```

Como puede apreciarse, en este caso se han ido ajustando esencialmente dos parámetros: el número de árboles empleados y la profundidad de estos, llegando a alcanzarse algo más de un 94% de acierto, ligeramente por encima del _random forest_. Para verificar la influencia de los dos parámetros citados no tneemos para que reprsentar gráficamente el modelo:

```{r}
plot(gbm)  # Explorar la influencia de parámetros en el modelo
```

## Importancia de las variables según RF

Un aspecto interesante de _random forest_ es que nos permite recuperar un ranking de importancia de las características (las variables), lo cual aporta cierto nivel de interpretabilidad a un modelo que está compuesto por cientos de árboles y que, en consecuencia, generaría una base de reglas inmesa. La función `varImpPlot()`, tomando como parámetro el mejor modelo producido anteriormente, nos ofrece este resultado:

```{r}
varImpPlot(rf$finalModel, 
           main = "Relevancia de características (RF)") 
```

Como podemos comprobar, la aparición en un correo electrónico de signos de exclamación, dólar y la palabra _remove_ son los que más influencia tienen a la hora de determinar si es legítimo o basura.

## Evaluación de los modelos con datos de test

Hasta ahora, en los apartados previos, hemos comprobado la estimación que el propio proceso de ajuste de parámetros ha llevado a cabo en cada _ensemble_. Dado que tenemos desde el inicio una partición con datos de test, veamos cuál es el rendimiento de cada modelo con datos que nunca ha visto antes:

### Bagging
```{r}
confusionMatrix(predict(bag, test), test$type)
```

### Random Forest
```{r}
confusionMatrix(predict(rf, test), test$type)
```

### Boosting
```{r}
confusionMatrix(predict(gbm, test), test$type)
```

De acuerdo a datos de test el _ensemble_ que mejor rendimiento alcanza es _random forest_, con más un de un 98% de acierto, seguido del _bagging_ y del _boosting_. Esto nos demuestra que la estimación llevada a cabo mediante validación interna, a pesar de usarse 10fcv y evaluar promedios de múltiples ejecuciones, no implica que el modelo con mejor estimación sea el que luego mejor se comporta con datos reales.

## Comparar los modelos entre sí

Podemos usar el método `resamples()` para llevar a cabo unan comparación, tanto numérica como gráfica, del rendimiento de los modelos basándonos en los datos de la estimación interna:

```{r}
diff <- resamples(list(Boosting = gbm, Bagging = bag, RF = rf))
summary(diff)
dotplot(diff)
```
Fijándonos en la columna `Max.` para la métrica `Accuracy`, en el resumen numérico, podemos observar que _boosting_ y _bagging_ apenas difieren en este caso, mientras que _random forest_ se comporta algo mejor. 

# Ejercicios adicionales

```{r}
mnist <- read.csv("datos/mnist_test.csv", header = FALSE)
dim(mnist)

sample(mnist$V1, 25)
```

## Examinamos el contenido de algunas muestras de MNIST
```{r}
prev.conf <- par(mfrow = c(5, 5), mai = c(0,0,0,0))
digits <- lapply(1:25, function(row) {
  digit <- matrix(as.numeric(mnist[row, 2:785]), ncol = 28, byrow = FALSE)
  image(digit[ , 28:1], col = gray(255:0 / 255), axes = FALSE)
})
par(prev.conf)
```

## Tomamos una pequeña parte de MNIST

```{r}
set.seed(4242)
datos <- mnist[sample(1:nrow(mnist), 1000), ]
datos[ , 1] <- as.factor(datos[ , 1])
levels(datos$V1) <- paste0("D", 0:9)
table(datos[ , 1]) # Número de muestras por clase
```

## Y la usamos para entrenar varios modelos

### kNN (_k Nearest Neightbors_)

```{r}
train5fcv <- trainControl(method = "cv", number = 5)

knn.mnist <- train(V1 ~ ., data = datos, method = "kknn", trControl = train5fcv)
knn.mnist
plot(knn.mnist)
```

### CART (_Classification and Regression Tree_)
```{r}
rpart.mnist <- train(V1 ~ ., data = datos, method = "rpart", trControl = train5fcv)
rpart.mnist
plot(rpart.mnist)
```

### SVM (_Support Vector Machine_)
```{r}
svm.mnist <- train(V1 ~ ., data = datos, method = "svmRadial", trControl = train5fcv)
svm.mnist
plot(svm.mnist)
```

### MLP (_MultiLayer Perceptron_)
```{r}

mlp.mnist <- train(V1 ~ ., data = datos, method = "mlp", trControl = train5fcv,
                   tuneGrid = expand.grid(size = c(30, 40, 50)))

mlp.mnist
plot(mlp.mnist)
```

### RF (_Random Forest_)
```{r}
rf.mnist <- train(V1 ~ ., data = datos, method = "rf", trControl = train5fcv)
rf.mnist
plot(rf.mnist)
```

### XGB (_Extreme Gradiente Boosting_)
```{r}
xgb.mnist <- train(V1 ~ ., data = datos, method = "xgbTree", trControl = train5fcv)
xgb.mnist
plot(xgb.mnist)
```


## Comparación de los modelos

```{r}
diff <- resamples(list(
  kNN = knn.mnist, CART = rpart.mnist, SVM = svm.mnist,
  MLP = mlp.mnist, RF = rf.mnist, XGB = xgb.mnist))

summary(diff)
dotplot(diff)
```


# Recursos de utilidad

- Libro [Análisis exploratorio y visualización de datos con R](https://fcharte.com/manuales/01_VisualizacionR/) del que soy autor y descargable desde mi  página web
- En la [sección Tutoriales](https://fcharte.com/tutoriales/) de mi página personal se ofrecen múltiples tutoriales sobre cómo realizar diversas tareas con R, por ejemplo:
   * [Redes neuronales con R](https://fcharte.com/tutoriales/00_annsconr/)
   * [Intercambio de datos entre Excel y R](https://fcharte.com/tutoriales/00_datos-excel-r/)
   * [Análisis de rendimiento de tareas en R](https://fcharte.com/tutoriales/00_rbenchmarking/)
- En mi [perfil de Github](https://github.com/fcharte?tab=repositories) hay distintos repositorios con contenidos de utilidad, desde paquetes R para acometer distintas tareas hasta el material que he creado para distintos cursos impartidos sobre R 
- Data Carpentry [Data Analysis and Visualization in R](https://datacarpentry.org/R-genomics/index.html), con ejercicios paso a paso sobre análisis de datos y visualización
- El libro [R for Data Science](https://r4ds.had.co.nz/) es una referencia gratuita (disponible en papel en Amazon) e imprescindible, escrito por el gurú de R Hadley Wickham
- En [MLR Playground](https://compstat-lmu.shinyapps.io/mlrPlayground/), una página web con R como motor, se pueden comparar las fronteras y resultados de distintos clasificadores aplicados a diversos conjuntos de datos