---
title: "Tutorial"
author: "Ricardo Aler"
date: "1 Agosto 2019"
output:
  pdf_document: default
  html_document: default
  word_document: default
---
# Tutorial sobre datos desbalanceados en mlR

## Introducción

Vamos a usar el conjunto de datos *ubIonosphere*, que tiene un ligero desbalanceo (proporción de positivos frente a negativos = `r 126/225`). Podemos ver que dado que la clase mayoritaria es un 64% de los datos, un clasificador trivial que predijera siempre "0" tendría un 64% de aciertos.

```{r}
data(ubIonosphere, package = "unbalanced")

# Comprobamos que la clase es un factor. Lo es, pero si no, habría que convertirla a factor.
class(ubIonosphere$Class)

(conteo <- table(ubIonosphere$Class))

conteo/sum(conteo)

#GRADO DE DESBALANCEO
conteo[2]/conteo[1]
```

El problema que tienen los datos desbalanceados es que las clases minoritarias se clasifican mal. Veamos si eso es así con k-vecinos. Tenemos que definir una tarea de clasificación. Es importante que digamos cual de las dos es la clase positiva (la minoritaria), que en este caso es la "1". También, cuando definamos el método para evaluar el modelo (con "Holdout" en este caso), tenemos que especificar que las particiones tienen que estar estratificadas. Por último, entrenamos / evaluamos el modelo con *resample*. Vamos a computar las medidas de error *tpr* y *tnr* (true positive rate y true negative rate), que son la tasa de aciertos de la clase positiva y negativa, respectivamente. También vamos a computar la tasa de aciertos (*acc*) y el "balanced accuracy", que es la media de *tpr* y *tnr*.

```{r}
library(mlr)

# Primero definimos la tarea. Nótese que definimos cual es la clase positiva
task_io <- makeClassifTask(data=ubIonosphere, target="Class", positive = "1")
learner_kknn <- makeLearner("classif.kknn")

# Al definir el método de evaluación, tenemos que especificar que las particiones sean estratificadas
metodo_evaluacion <- makeResampleDesc("Holdout", stratify=TRUE) 
set.seed(0)
particion_datos <- makeResampleInstance(metodo_evaluacion, task_io)
                         
set.seed(0) 
errores_kknn <- resample(learner_kknn, task_io, particion_datos, measures = list(acc, bac, tpr, tnr))
```

Vemos que la clase negativa (mayoritaria) se aprende mucho mejor que la positiva: tnr = `r errores_kknn$aggr["tnr.test.mean"]` y tpr = `r errores_kknn$aggr["tpr.test.mean"]`. Aún así, dado que *acc* (83%) supera el 65% de aciertos, el clasificador construido con *knn* es mejor que un clasificador trivial. Vemos también que el error *bac* representa mejor que *acc* la precisión del clasificador, sobre todo con respecto a la clase positiva (minoritaria).

```{r}
errores_kknn$aggr
```
La matriz de confusión muestra similares resultados 

```{r}
calculateConfusionMatrix(errores_kknn$pred, relative = TRUE)
```
## Mostrando curvas ROC

Vamos a hacer una curva ROC para *kknn*. Para ello es necesario hacer predicciones probabilísticas. Aquí las podemos ver:

```{r}

learner_kknn_prob <- makeLearner("classif.kknn", predict.type = "prob")
model_kknn_prob <- train(learner_kknn_prob, task_io, subset = particion_datos$train.inds[[1]])
predicciones <- predict(model_kknn_prob,  task_io, subset = particion_datos$test.inds[[1]])
head(predicciones$data)

```
Pero en este momento, no nos interesan las predicciones, vamos a calcular los errores directamente con *resample*. Ahora que hacemos predicciones probabilísticas, podemos aprovechar para calcular también el área bajo la curva (auc), como se puede ver abajo.

```{r}
set.seed(0) 
errores_kknn <- resample(learner_kknn_prob, task_io, particion_datos, measures = list(acc, bac, tpr, tnr, auc))
```
Para generar la curva ROC, primero tenemos que generar cierta información con la función *generateThreshVsPerfData*, a partir de los errores computados anteriormente. Necesitaremos *fpr* y *tpr* para construir la curva ROC, pero podemos generar información también para otras medidas de error, como *tnr* y *bac*. Vemos que la información generada no es mas que un data.frame que contiene una línea para cada valor del threshold.

```{r}
df <- generateThreshVsPerfData(errores_kknn, measures = list(fpr, tpr, tnr, bac))
head(df$data)
```
Ahora imprimimos la curva ROC. Normalmente, imprime sólo la línea, pero si queremos que dibuje también los puntos, usamos *geom_point()*. mlR genera los plots mediante *ggplot2*.

```{r}
library(ggplot2)
plotROCCurves(df, measures=list(fpr, tpr)) + geom_point()
```
Se puede usar *plotly* para hacer plots ligeramente interactivos. Podemos ver que podríamos mejorar bastante el *tpr*, sin empeorar demasiado el *tnr*. Originalmente obteníamos un *tpr* de 0.59, con un *tnr* de 0.96 (o lo que es lo mismo, un fpr = 1-tnr = 0.04). Vemos que cambiando el threshold, podríamos subir el *tpr* a 0.88 y el *tnr=1-fpr* sólo bajaría a 1-0.05 = 0.95.

```{r}
library(ggplot2)
library(plotly)
gg <- plotROCCurves(df, measures=list(fpr, tpr)) + geom_point()

ggplotly(gg)
```

## Ajustando el threshold a mano

Desgraciadamente, en la curva ROC, no se puede ver el threshold que tendríamos que utilizar para obtener dicho resultado. Vamos a generar otra curva, en cuyo eje x viene el threshold, y en el eje y vienen los distintos errores. Haremos el plot con *ggplot2*. Para ello, primero tenemos que poner los datos en modo *largo* (long)

```{r}
# Primero p
df_long <- tidyr::gather(df$data, variable, value, -threshold)
head(df_long)
```

Ahora hacemos el plot. Vamos a hacerlo interactivo. Podemos ver que para un threshold de 0.01 se alcanza un buen equilibrio entre *tpr* y *tnr*. Podemos usar para medirlo la media de los dos *bac*.

```{r}
library(ggplot2)
df_long_tpr_tnr_bac <- df_long[df_long$variable %in% c("tpr", "tnr", "bac"), ]
gg <- ggplot(df_long_tpr_tnr_bac, aes(x=threshold, y=value, col=variable )) + geom_point() + geom_line()
ggplotly(gg)

```
Podemos calcular directamente el threshold, como aquel lugar donde se alcanza el máximo *bac*, así:

```{r}
th <- df$data[which.max(df$data$bac), "threshold"]
th
```

Vamos a fijar dicho threshold con *setThreshold* al hacer las predicciones.

```{r}
preds_kknn_th <- setThreshold(errores_kknn$pred, th)
performance(preds_kknn_th, measures=list(tpr, tnr, bac))
```
## Ajustando el threshold automáticamente

También podemos ajustarlo automáticamente con *tuneThreshold* Hay que fijarse que estamos encontrando el threshold que optimiza *bac*

```{r}
# El threshold óptimo es
(th <- tuneThreshold(errores_kknn$pred, measure=list(bac), nsub=100))

preds_kknn_th <- setThreshold(errores_kknn$pred, th$th)
errores_kknn_tuned <- performance(preds_kknn_th, measures=list(tpr, tnr, bac))
```

Vemos que obtenemos el mismo *bac* que cuando ajustamos el threshold a mano.

```{r}
errores_kknn_tuned
```
## Ajustando el threshold como cualquier otro hiper-parámetro

En los ajustes de threshold que hemos hecho anteriormente había un aspecto que no era del todo correcto. Tanto las curvas ROC que visualizamos, como el ajuste visual o automático del threshold, lo hicimos con las predicciones en la partición de test, lo cual no es correcto, dado que estamos usando el conjunto de test para ajustar un parámetro del modelo (en este caso, el threshold). Viene bien para visualizar como varían los distintos errores (tpr, tnr, bac, etc.) a medida que cambia el threshold. Pero el ajuste del threshold habría que hacerlo como si fuera cualquier otro hiper-parámetro (la manera más sencilla es con un conjunto de validación, distinto al de test). Eso es lo que vamos a hacer a continuación. En mlR, esto se puede hacer usando el ajuste de hiper-parámetros que ya conocemos, excepto que en *makeTuneControlGrid*, tenemos que añadir los argumentos *tune.threshold* y *tune.threshold.args*, como vemos en el siguiente código. Nötese que estamos pidiendo que ajuste el threshold optimizando *bac* (podríamos optimizar cualquier otra medida).

Si quisieramos, podríamos ajustar a la vez el hiper-parámetro *k*, pero en este caso, sólo ajustaremos el threshold. Esa es al razón por la que hemos fijado la lista de valores de *k* a un único valor, el 7 (pero podríamos poner una lista de valores, por ejemplo, de 1 a 10). 

```{r}
getParamSet(learner_kknn_prob)

ps = makeParamSet(
  makeDiscreteParam("k", values = 7)
)

control_grid <- makeTuneControlGrid(tune.threshold = TRUE, tune.threshold.args = list(measure=list(bac)))
evaluacion_grid <- makeResampleDesc("Holdout", stratify = TRUE)

learner_ajuste_kknn <- makeTuneWrapper(learner_kknn_prob, resampling = evaluacion_grid, par.set = ps, control = control_grid, measures = list(bac))

set.seed(0)
errores_ajuste_kknn <- resample(learner_ajuste_kknn, task_io, particion_datos, measures = list(acc, bac, tpr, tnr), extract = getTuneResult, models = TRUE)
```
Vemos que el threshold es ligeramente distinto al de antes, y también los errores. Esta diferencia es debida a que antes optimizabamos el threshold directamente con el conjunto de test (que ya sabemos que no es correcto), y ahora lo estamos optimizando con un conjunto de validación. Una vez optimizado el threshold, evaluamos el modelo con este threshold sobre el conjunto de test.

```{r}
errores_ajuste_kknn$models[[1]]$learner.model$opt.result$threshold
errores_ajuste_kknn
```
Podemos aprovechar para ajustar también otros hiper-parámetros (como k), además del threshold:

```{r}

ps = makeParamSet(
  makeDiscreteParam("k", values = 2:10)
)

control_grid <- makeTuneControlGrid(tune.threshold = TRUE, tune.threshold.args = list(measure=list(bac)))
evaluacion_grid <- makeResampleDesc("Holdout", stratify = TRUE)

learner_ajuste_kknn <- makeTuneWrapper(learner_kknn_prob, resampling = evaluacion_grid, par.set = ps, control = control_grid, measures = list(bac))

set.seed(0)
errores_ajuste_kknn <- resample(learner_ajuste_kknn, task_io, particion_datos, measures = list(acc, bac, tpr, tnr), extract = getTuneResult, models = TRUE)
```



```{r}
errores_ajuste_kknn$models[[1]]$learner.model$opt.result$threshold
errores_ajuste_kknn$models[[1]]$learner.model$opt.result$x
errores_ajuste_kknn
```
## Usando SMOTE para abordar problemas de muestras desbalanceadas

En primer lugar, vamos a rebalancear la muestra "a mano", haciendo que las dos clases tengan el mismo número de instancias. Si la clase mayoritaria tuviera el doble de instancias que la minoritaria, habría que duplicar el número de instancias positivas (o sea, un *rate* de 2). En este caso, usaremos un *rate* de instancias mayoritaria / instancias minoritaria. Para usar SMOTE, montaremos un método híbrido, cuya primera fase hace el rebalanceo con SMOTE, y cuya segunda fase aplica *knn* sobre la muestra desbalanceada.

```{r}
tabla <- table(getTaskTargets(task_io))
tabla

rate <- tabla["0"]/tabla["1"]
rate

#task_m_smote <- smote(task_m, rate = rate)
#table(getTaskTargets(task_m_smote))

learner_smote_kknn = makeSMOTEWrapper(learner_kknn_prob, sw.rate = rate)

set.seed(0) 
errores_smote_knn <- resample(learner_smote_kknn, task_io, particion_datos, measures = list(acc, bac, tpr, tnr), models = TRUE)

```

Vemos que con este *rate*, conseguimos un error parecido al que obtuvimos antes ajustando el threshold.

```{r}
errores_smote_knn$aggr
```
En lugar de poner el *rate* a mano, podemos ajustarlo automáticamente, dado que es un hiper-parámetro del método híbrido *learner_smote_kknn*. Usaremos una validación cruzada de 3 folds para ajustar *rate*.

```{r}
ps = makeParamSet(
  makeDiscreteParam("sw.rate", values = seq(1, 10, 1))
)

evaluacion_grid <- makeResampleDesc("CV", iters=3, stratify = TRUE)
control_grid <- makeTuneControlGrid()

learner_smote_kknn = makeSMOTEWrapper(learner_kknn_prob)

learner_ajuste_smote_kknn <- makeTuneWrapper(learner_smote_kknn, resampling = evaluacion_grid, par.set = ps, control = control_grid, 
                                              measures = list(bac, tpr, acc, tnr), show.info = FALSE)

set.seed(0) 
errores_ajuste_smote_kknn <- resample(learner_ajuste_smote_kknn, task_io, particion_datos, measures = list(acc, bac, tpr, tnr), extract = getTuneResult, models=TRUE)
```
Vemos que conseguimos un *bac* superior al que hemos ido obteniendo hasta el momento. Vemos que para obtener el mejor resultado, fue necesario utilizar una replicación cuadruple de la clase minoritaria.

```{r}
errores_ajuste_smote_kknn

errores_ajuste_smote_kknn$models[[1]]$learner.model$opt.result

  

```
Puede ser interesante ver como cambia el error al cambiar el *rate*. Vemos que multiplicar por 8 las instancias de la clase minoritaria tiene un efecto beneficioso.

```{r}
efectos <- generateHyperParsEffectData(errores_ajuste_smote_kknn)
head(efectos)
plotHyperParsEffect(efectos, x="sw.rate", y="bac.test.mean", loess.smooth = TRUE)

```


Además de ajustar el *rate* de SMOTE, también podemos optimizar el threshold (esto lo podemos hacer al ajustar hiper-parámetros, siempre que el método permita predecir probabilidades, que es el caso de *knn*).

```{r}
ps = makeParamSet(
  makeDiscreteParam("sw.rate", values = seq(1, 10, 1))
)

evaluacion_grid <- makeResampleDesc("CV", iters=3, stratify = TRUE)
control_grid <- makeTuneControlGrid(tune.threshold = TRUE, tune.threshold.args = list(measure=list(bac)))

learner_smote_kknn = makeSMOTEWrapper(learner_kknn_prob)

learner_ajuste_smote_kknn <- makeTuneWrapper(learner_smote_kknn, resampling = evaluacion_grid, par.set = ps, control = control_grid, 
                                              measures = list(bac, tpr, acc, tnr), show.info = FALSE)

set.seed(0) 
errores_ajuste_smote_kknn <- resample(learner_ajuste_smote_kknn, task_io, particion_datos, measures = list(acc, bac, tpr, tnr), extract = getTuneResult, models=TRUE)
```
Parece que usando un *rate* de 10 y un threshold de 0.35, podemos mejorar un poco más el *bac*.

```{r}
errores_ajuste_smote_kknn

errores_ajuste_smote_kknn$models[[1]]$learner.model$opt.result

```
```{r}
efectos <- generateHyperParsEffectData(errores_ajuste_smote_kknn)
head(efectos)
plotHyperParsEffect(efectos, x="sw.rate", y="bac.test.mean", loess.smooth = TRUE)

```

Incluso podemos ajustar el mismo tiempo el hiper-parámetro *k*, aunque esto va a tardar más tiempo.

```{r}
library("parallelMap")
parallelStartSocket(4)

ps = makeParamSet(
  makeDiscreteParam("sw.rate", values = seq(1, 10, 1)),
  makeDiscreteParam("k", values = seq(1, 10, by=2))
)

evaluacion_grid <- makeResampleDesc("CV", iters=3, stratify = TRUE)
control_grid <- makeTuneControlGrid(tune.threshold = TRUE, tune.threshold.args = list(measure=list(bac)))

learner_smote_kknn = makeSMOTEWrapper(learner_kknn_prob)

learner_ajuste_smote_kknn <- makeTuneWrapper(learner_smote_kknn, resampling = evaluacion_grid, par.set = ps, control = control_grid, 
                                              measures = list(bac, tpr, acc, tnr), show.info = TRUE)

set.seed(0) 
errores_ajuste_smote_kknn <- resample(learner_ajuste_smote_kknn, task_io, particion_datos, measures = list(acc, bac, tpr, tnr), extract = getTuneResult, models=TRUE)

parallelStop()
```
Pero en este caso, ajustar *k* no permite mejorar *bac*.

```{r}
errores_ajuste_smote_kknn

errores_ajuste_smote_kknn$models[[1]]$learner.model$opt.result

```
```{r}
efectos <- generateHyperParsEffectData(errores_ajuste_smote_kknn)
head(efectos)
plotHyperParsEffect(efectos, x="sw.rate", y="bac.test.mean", loess.smooth = TRUE)
plotHyperParsEffect(efectos, x="k", y="bac.test.mean", loess.smooth = TRUE)
plotHyperParsEffect(efectos, x="k", y="sw.rate", z="bac.test.mean", plot.type = "heatmap")

```

Nota: hay métodos que permiten dar pesos a las instancias. En ese caso, una aproximación más rápida que SMOTE podría ser usar *makeWeightedClassesWrapper* en lugar de *makeSMOTEWrapper*. En este caso, el nombre del hiper-parámetro que mide el peso (o importancia) de los datos de la clase minoritaria frente a los de la mayoritaria es *wcw.weight*. Desgraciadamente, **kknn** no permite usar pesos con las instancias (pero *rpart* si).