68_SDM_Algorithms.Rmd

# Algoritmos de modelos de distribucion de especies

Introducción a varios algoritmos disponibles para modelar distribución de especies.  
  
[Presentation (PDF)](./course_material/C6_Different algorithms.pdf)  
  
## Ejercicios - Práctica con modelos

Para este ejercicio vamos a utilizar un set de datos de superhéroes. De seguro que se sorprenderán cuando comiencen a explorar los datos ambientales.  
  
El objetivo es correr modelos que nos permitan responder dos preguntas:  
  
1) CLASIFICACIÓN: ¿Existe un sesgo hacia la armonización? Es decir, ¿somos capaces de predecir si un personaje es bueno o malo?  
2) REGRESIÓN: ¿Son los hombres de los cómics más fuertes que las mujeres? Es decir, ¿somos capaces de predecir la fuerza total del personaje? Además, ¿es el género uno de los predictores más importantes o hay alguna otra variable que es mejor para predecir la fuerza total?  
  
Aquí haremos lo siguiente:  
  
1) Explorar y limpiar datos  
2) Usar el método de codificación "one-hot" para incluir muchas categorías en el modelo  
3) Construir un árbol de decisión y visualizar resultados. Reconstruir este árbol pero con configuración diferente  
4) Correr un modelo con `gbm.step` y visualizar las figuras de entrenamiento y prueba, así como crear una figura de la importancia de las variables usadas en el modelo  
5) Correr el mismo modelo pero con random forests y contrastar las diferencias en los resultados      
  
***Escoge un problema de clasificación o regresión. La solución para el problema de clasificación está incluida en este documento por si necesitas ayuda. Aquí vamos a correr el modelo GBM solamente para el problema de regresión, pero también puedes probarlo con random forest (RF)***  
  
Empecemos llamando a las librerías que necesitaremos  

```{r sdmalg1, warning=FALSE}
library(tidyverse)
## Esta librería incluye el comando gbm.step
library(dismo)
## Esta librería incluye randomForest
library(randomForest)  
## Esta librería incluye una serie de funciones útiles para los árboles de 
#decisión
library(rpart)   
library(rpart.plot)
library(gbm)
library(GGally)
```
  
### Cargando y limpiando datos

Una vez llamadas las librerías relevantes, vamos a subir el set de datos de superhéroes y los limpiaremos antes de poder iniciar con los modelos.  
  
```{r sdmalg2, echo = FALSE}
## Aqui leemos los datos de superheroes
Supes <- read_csv('docs/course_material/exercises/data/SuperheroDataset.csv')
```
  
Para este ejercicios vamos a utilizar las siguientes columnas: Name, Intelligence, Strength, Speed, Durability, Power, Combat, Gender, Race, Creator, Alignment, Total Power. Usamos `tidyverse` para extraer estas columnas.   
  
```{r sdmalg4}
SuperData <- Supes |>
  dplyr::select(Name, Intelligence:Combat, Alignment:Race, Creator, TotalPower)
```
  
#### Explorando los datos
  
Si exploramos los datos podemos ver algunos elementos interesantes de este set de datos. Usamos `ggplot2` (o `plot` de `R` base si prefieres). 

```{r sdmalg5, warning=FALSE,message=FALSE}

SuperData |>
  filter(!is.na(Intelligence))|> 
  ggplot(aes(x = Intelligence, y = Strength)) + 
  geom_point() + 
  geom_smooth(method = lm) +
  theme_bw()

## Mirando a la linea parece que existe un valor extremadamente bajo de 
# inteligencia. Probablemete puntos atipicos. Veremos a quien le pertenecen 
# estos valores, y que podria ayudar a determinar si deberiamos incluirlos en 
# nuestro analisis

SuperData |>
  filter(!is.na(Combat))|> 
  ggplot(aes(x = Combat, y = Strength)) + 
  geom_point() + 
  geom_smooth(method = lm) + 
  theme_bw()

## Podemos tambien filtrar datos simplemente para explorar el set de datos 
# usando condiciones especificas, como por ejemplo:
SuperData |>
  filter(!is.na(Combat), Gender == 'Female')|> 
  ggplot(aes(x = Combat, y = Strength)) + 
  geom_point() + 
  geom_smooth(method = lm) + 
  theme_bw()

## Esta es una buena librería para visualizar datos rapidamente en lugar de 
# crear gráficos unicos con codigo

quant_df <- SuperData |> 
  dplyr::select(Intelligence:Gender, TotalPower)
ggpairs(quant_df, progress = FALSE)

```
  
Nota que `total power` está altamente correlacionado con `strength` (fuerza), `speed` (velocidad), `durability` (durabilidad) y `power` (poder). Esto podría influenciar cómo interpretamos el problema de regresión, ya que la correlación de estas variables pueden interferir con la señal de influencia de las variables exploradas. Esto puede interferir en nuestra habilidad de responder nuestras preguntas.  

Si has explorado el set de datos quizás hayas notado lo siguiente:  
  
1)  `Race` (raza) and `Creator` (creador) tienen muchas categorías.  
2)  Existen muchos creadores que solo fueron mencionados una vez, y quizás no sea apropriado incluirlos todos.  
  
```{r sdmalg6}
## Cantidad de categorias unicas en raza
length(unique(SuperData$Race))
## Lo mismo pero para creador
length(unique(SuperData$Creator))  

## Explremos a Creator:
unique(SuperData$Creator)

SuperData |> 
  count(Creator)

```

Esta columna tiene muchas categorías y probablemente no sea apropriado incluir todas en el análisis. Por ejemplo, Ian Fleming (creador de James Bond) está asociado a un solo personaje. Este set de datos también tiene personajes de Star Trek y muchas otras películas o series de televisión. Quizás podamos solo incluir personajes de cómics como Marvel Comics, DC Comics, Dark Horse Comics, e Image Comics.   
  
***Nota el uso de `%in%`***

```{r sdmalg7}
FilteredSupes <- SuperData |>
  dplyr::filter(Creator %in% 
                  c('Marvel Comics', 'DC Comics', 'Dark Horse Comics', 
                    'Image Comics')) |>
  dplyr::filter(!is.na(Intelligence))

## la ultima linea remueve NA del set de datos. Este paso no es estrictamente 
# necesario con arboles de decision

FilteredSupes

## ¿Cuantas razas tenemos?
length(unique(FilteredSupes$Race))
```
  
Exploremos los datos otra vez.  
  
```{r sdmalg8, warning = FALSE, message = FALSE}
quant_df <- FilteredSupes |> 
  dplyr::select(Intelligence:Gender, TotalPower)
ggpairs(quant_df, progress = FALSE)
```
  
No mucho ha cambiado porque no removimos muchas filas. Esto es bueno porque mantenemos la mayor parte de nuestros datos y esperemos que las conclusiones finales no sean muy diferentes. Pero el beneficio de remover esas pocas filas es que nos ayuda a simplificar el proceso de limpieza de nuestros datos para este ejercicio. Sin embargo, aún tenemos muchas categorías bajo `Race`.  
  
```{r sdmalg9}
FilteredSupes |> 
  count(Race)
```
  
Mirando en más detalle a `Race`, vemos varias categorías con una sola observación. Esto puedo ocasionar problemas cuando hagamos validación cruzada, así que es preferible removerlos.  
  
```{r sdmalg10}
FilteredSupes <- FilteredSupes |> 
  add_count(Race) |> 
  filter(n > 1)

FilteredSupes |>
  count(Race)
```

### Codificación One-hot

`Race` tiene muchas categoría y esto podrías ocasionar que se pierda la señal de las variables utilizadas como predictores por los árboles de decisión.  Para evitar problemas vamos a hacer uso de la codificación one-hot, la cual nos permitirá mantener el mayor número de categorías posibles. La codificación one-hot nos permite crear categorías utilizando variables binarias, y esto se puede hacer de manera sencilla con el `tidyverse`.  
  
```{r sdmalg11}
OOEdata <- FilteredSupes |> 
        separate_rows(Race)|> 
        mutate(count = 1) |> 
        spread(Race, count, fill = 0, sep = "_")

OOEdata
```
  
Ahora podemos explorar las categorías fácilmente.

***Recuerda que necesitamos convertir estas categorías a factores. De otra manera serán considerados como números enteros***  
 
La librería ***magrittr*** nos puede ayudar a hacer esto rápidamente.  

```{r sdmalg12, message = FALSE, warning = FALSE}
library(magrittr)

columns <- names(OOEdata)[12:length(names(OOEdata))]
OOEdata <- OOEdata %<>% mutate_at(columns, funs(factor(.)))


quant_df <- OOEdata |>
  dplyr::select(Intelligence:TotalPower, Race_Human, Race_Cyborg, Race_Mutant)

ggpairs(quant_df, progress = FALSE)

```

Este gráfico no se muestra muy bien en el `R Markdown`, pero si excluimos algunas variables podríamos obtener un mejor gráfico. También es posible usar el visualizador de gráficos de R Studio, el cual tiene un botón de zoom que nos puede ayudar a explorar el gráfico de mejor manera.   
  
### Los modelos  

#### CART
  
Comencemos construyendo un solo árcbol. Puedes usar cualquier variable aquí.  

##### Ejemplo de clasificación  
  
```{r sdmalg13}
## Vamos a remover los personajes 'neutrales' de nuestro set de datos y vamos
# a intentar predicir si los personajes son buenos o malos. Tambien vamos a 
# seleccionar solo unas pocas columnas para simplificar nuestro trabajo.

modeldat <- OOEdata |> 
  dplyr::filter(Alignment!='neutral') |>
  dplyr::select(Intelligence:TotalPower, Race_Animal, Race_Human, Race_Demon,
                Race_God, Race_Asgardian, Race_Mutant, Race_Android)

form <- as.formula(paste("Alignment ~ ", 
                         paste(names(modeldat), collapse = '+')))

form <- update(form, . ~ . -Alignment) 

## Revisa help(rpart.control) para obtener una lista de parametros que puedes 
# cambiar en el modelo. Luego intentar cambiar los parametros para que te des 
# una idea del impacto que tienen al crear un solo arbol

## Para correr un modelo de regresion, cambia "method" a 'anova'
treemod <- rpart(form, method = 'class', data = modeldat, cp = 0.001, 
                 maxdepth = 6)

rpart.plot(treemod, cex = 0.7)

```
  
#### GBM.step
  
Nos saltaremos un paso aquí y usaremos el algoritmo `gbm.step`, el cual es una extensión de `GBM` (Generalized Boosted regression Modelling también conocido como boosted regression trees). `gbm.step` calculará el árbol óptimo que deberemos usar para obtener un modelo mejor poder de predicción. Para esto utilizará validación cruzada, en vez de contruir varios árboles y dejar que el usuario decida por si mismo.  
  
Utilizaremos un modelo de clasificación como el de más arriba. Pero debemos tomar en cuenta que `gbm.step` necesita el índice las columnas en vez de sus nombres.
  
Notas adicionales de `gbm.step`:
***`gbm.step` no está diseñado para considerar a múltiples categorías meta (variables dependientes)***
***`gbm.step` acepta datos solo en un `data frame`, no es posible usar un `tibble`. Además la categoría a ser predicha debe ser de tipo `factor`***
***`gbm.step` también necesita una categoría meta forma binaria (0/1)***

```{r sdmalg14}
modeldf <- modeldat |> 
  mutate(Alignment_targ = case_when(Alignment == "good" ~ 0, 
                                    Alignment == "bad" ~ 1 )) |>
  mutate(across(c("Alignment_targ", "Gender", "Creator"), factor)) |> 
  as.data.frame()

modeldf$Alignment_targ <- as.numeric(as.character(modeldf$Alignment_targ))

gbmModel <- gbm.step(data = modeldf, gbm. = c(1:6,8:17), gbm.y = 18, 
                     family = "bernoulli", tree.complexity = 10, 
                     learning.rate = 0.001, bag.fraction = 0.8, 
                     max.trees = 8000)
```
  
Para obtener un resumen de los resultados que muestren una clasificación de la importancia de las variables usa la línea de abajo:
  
```{r sdmalg15}
summary(gbmModel)
```
  
Puedes también crear gráficos de dependencia parcial, y además puedes explorar las interacciones entre variables con las líneas de código de abajo  
  
```{r sdmalg16}
gbm.plot(gbmModel, n.plots = 5, write.title = F)

## el gbm.perspec produce un gráfico de inteligencia y velocidad contra los 
# valores ajustados (z)
gbm.perspec(gbmModel, 1, 3, z.range = c(0.15, 0.6), theta = 205)

```
  
#### Random forest
  
De muchas maneras, un random forest es mucho más flexible que `gbm.step`. Por ejemplo, no estamos forzados a utilizar índices en vez de nombres de columnas.  
  
Este modelo puede ser configurado de manera muy similar al modelo `CART` que vimos antes.   
  
```{r sdmalg17}
modeldat <- OOEdata |> 
  dplyr::filter(Alignment!="neutral") |>
  dplyr::select(Intelligence:TotalPower, Race_Animal, Race_Human, Race_Demon, 
                Race_God, Race_Asgardian, Race_Mutant, Race_Android)

## Asegurate que las categorias esten en formato factor
modeldat <- modeldat |> 
  mutate(Alignment_targ = case_when(Alignment == 'good' ~ 0,
                                    Alignment == 'bad' ~ 1 )) |>
  mutate(across(c("Alignment_targ", "Gender", "Creator"), factor))

form <- as.formula(paste("Alignment_targ ~ ", paste(names(modeldat), 
                                                    collapse = '+')))
form <- update(form, . ~ . -Alignment) 
form <- update(form, . ~ . -Alignment_targ) 

rf.model <- randomForest(data = data.frame(modeldat), form, mtry = 3, 
                         maxnodes = 10)
importance(rf.model)
partialPlot(rf.model, data.frame(modeldat), Intelligence)
partialPlot(rf.model, data.frame(modeldat), Gender)
```
  
***Los gráficos de dependencia parcial se enfoca en la PRIMERA clase. Esto quiere decir que si usamos clases 0 y 1, los gráficos deberán ser interpretados en referencia al grupo 0. Pero podemos cambiar esto con el argumento `which.class` en la función `partialPlot`***
  
### Problemas de regresión
  
Ahora vamos a utilizar `gbm.step` para resolver un problema de regresión. Si quieres usar random forest, también es posible y te alentamos a intentarlo. Puedes pedir ayuda a los instructores o si prefieres puedes usar Google, ya que hay muchos recursos disponibles en línea.  
  
```{r sdmalg18}
modeldat <- OOEdata |> 
  dplyr::filter(Alignment != "neutral") |>
  dplyr::select(Intelligence:TotalPower, Race_Animal, Race_Human, Race_Demon, 
                Race_God, Race_Asgardian, Race_Mutant, Race_Android)

modeldf <- modeldat |> 
  mutate(across(c("Alignment", "Gender", "Creator"), factor)) |> 
  as.data.frame()

## usa names(modeldf) para encontrar los indices

gbmModel <- gbm.step(data = modeldf, gbm.x = c(1,6:9,11:17), gbm.y = 10,
                     family = "laplace", tree.complexity = 20,
                     learning.rate = 0.001, bag.fraction = 0.5, 
                     max.trees = 8000)

summary(gbmModel)

gbm.plot(gbmModel, n.plots = 5, write.title = F)

## En este caso estamos interesados en el genero, asi que vamos a graficar esto.

gbm.plot(gbmModel, variable.no = 4, plot.layout = c(1, 1))
```
  
Podemos ver que la fuerza de un personaje está sesgada de acuerdo al género, con más héroes que heroínas con mucho poder. Sin embargo, vale recalcar que el poder predictivo de nuestro modelo no es muy alto (basado en el resumen estadístico).