PracticaAC.Rmd

---
title: "Análisis de Datos y Machine Learning con R (Caret)"
author: ["Eva Pina Dubovtseva (G1.4)", "Jesús Sánchez Pardo(G1.5)", "Jiahui Lin(G1.2)"]
output: 
  pdf_document:
    toc: TRUE
    toc_depth: 5
    number_sections: true
date: "10-12-2024"
editor_options: 
  markdown: 
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

#Cargar el espacio de trabajo
load("ProyectoAC.Rdata")
```
\newpage
# Librerías empleadas

```{r librerías, message=FALSE}
library(lattice)
library(GGally)
library(ggplot2)
library(reshape2)
library(gridExtra)
library(caret)
library(mlbench)
library(rpart)
library(rpart.plot)
library(nnet)
library(randomForest)
library(gbm)
library(NeuralNetTools)
library(tidyr)
library(dplyr)
require(plotrix)
```

# Análisis Explotario de Datos

## Carga del dataset

Lo primero de todo, vamos a cargar en la variable "credit" todo el dataset.
```{r carga-datos}
# Cargar dataset
credit <- read.table("crx.data", sep=",", header=FALSE, na.strings="?")
```

## División en Train y Test

Para poder entrenar los modelos, es **esencial** **dividir los datos** que 
tenemos **en un conjunto de entrenamiento**, **y otro de test** 
para evaluarlo en el modelo final.

Se van a utilizar unos índices preestablecidos para ello.
```{r division}
# Dividimos el dataset en train y test
credit.trainIdx <- readRDS("credit.trainIdx.rds")
credit.Datos.Train <- credit[credit.trainIdx, ]
credit.Datos.Test <- credit[-credit.trainIdx, ]
```

Vamos a visualizar el dataset entero:

```{r visualizacion-inicial}
# Estructura del dataset
head(credit)

```

Como podemos ver, consta de 16 columnas. **Al estar codificado, es complicado**
**interpretar qué es cada variable**. No obstante, claramente la última columna 
(**V16**) está relacionado con si el crédito ha sido aprobado o no. Está última,
**será nuesta columna objetivo a lo largo de toda la práctica**.

Con el comando summary() vamos a tener un resumen estadístico de cada variable. 
De esta manera, podremos tener una visión superficial de las características de 
los datos con los que estamos trabajando. También lo haremos con 
los conjuntos de train y test.

```{r}
summary(credit)
summary(credit.Datos.Train)
summary(credit.Datos.Test)
```


La salida nos muestra que las columnas **V2**, **V3**, **V8**, **V11** y **V15** son numéricas. Para ellas, se han calculado valores como la **media**, la **mediana**, entre otros. Si analizamos bien cada una, podemos destacar:

-   **V2**: Tiene valores que varían entre **13.75** y **80.25**. La **mediana** de **28.46** sugiere que los valores están **sesgados hacia la izquierda** (menores valores).

-   **V3**: Presenta una **mediana de 2.75** y un máximo de **28.00**, lo que indica que la mayoría de los valores están en un **rango bajo**.

-   **V15**: Tiene valores muy amplios, desde **0** hasta **1,000,000.0**, con una **media de 1017.4**. Esto puede indicar la presencia de **outliers**.

Por otro lado, también se observa que **en el conjunto de train, y todo el dataset, hay valores nulos en algunas columnas**. En la siguiente sección, se tratarán estos nulos para poder visualizar y analizar los datos completos, y garantizar que los modelos sean entrenados adecuadamente con el conjunto de train. **El conjunto de test no presenta valores nulos**, por lo que no necesita tratamiento.

Finalmente, se procede a **guardar las columnas numéricas y categóricas**. Utilizando el comando **sapply()**, se aplica la función **is.numeric** o **is.factor** para identificar los tipos de columnas y guardar los resultados en vectores correspondientes.


```{r}
# Identificar tipos de columnas
numeric_columns <- sapply(credit, is.numeric)
categorical_columns <- sapply(credit, is.factor)

# Guardaremos también los nombres de las columnas
numeric_columns2 <- names(credit)[sapply(credit, is.numeric)]
categorical_columns2 <- names(credit)[sapply(credit, function(x) is.character(x) || is.factor(x))]

```

## Tratamiento de valores nulos

Como hemos visto con el summary, tenemos datos nulos. Vamos a tratarlos
primero con tal de tener los datos actualizados a nuestro análisis.

Con el siguiente comando, podemos ver cuántos datos nulos hay en cada
variable.

```{r}
sapply(credit, function(x) sum(is.na(x)))

sapply(credit.Datos.Train, function(x) sum(is.na(x)))
```

Al ser un dataset **no muy grande**, vamos a optar por **rellenar los NAs** con la **mediana** para las variables de valores **numéricos**, y el **más frecuente** si es un **factor**. Además, de que **no son demasiados los valores faltantes**, por lo que el **relleno de los NAs no debería introducir mucho sesgo**.


```{r categorias, fig.width=5, fig.height=3}
# Convertir columnas categóricas a factor
credit[categorical_columns2] <- lapply(credit[categorical_columns2], as.factor)
credit.Datos.Train[categorical_columns2] <- lapply(credit.Datos.Train[categorical_columns2], as.factor)
credit.Datos.Test[categorical_columns2]  <- lapply(credit.Datos.Test[categorical_columns2], as.factor)

```

```{r valores-faltantes}
# Número total y por columna de valores faltantes
sum(is.na(credit))

# Rellenar valores faltantes
fixNA_credit <- function(data.toFix) {
  for (i in seq_along(data.toFix)) {
    if (is.factor(data.toFix[[i]])) {
      ttx <- table(data.toFix[[i]], useNA = "no")
      data.toFix[[i]][is.na(data.toFix[[i]])] <- names(ttx)[which.max(ttx)]
    } else if (is.numeric(data.toFix[[i]])) {
      data.toFix[[i]][is.na(data.toFix[[i]])] <- median(data.toFix[[i]], na.rm = TRUE)
    }
  }
  return(data.toFix)
}

# Aplicamos la fórmula en credit y en train
credit <- fixNA_credit(credit)
credit.Datos.Train <- fixNA_credit(credit.Datos.Train)

# Verificar nuevamente valores faltantes
sapply(credit, function(x) sum(is.na(x)))
sapply(credit.Datos.Train, function(x) sum(is.na(x)))
```

Como se puede ver, **ya no quedan nulos en el dataset**, **ni en train**.

## Análisis univariable

Al no conocer nada sobre el conjunto, hemos decidido hacer un análisis de todas las variables para ver la distribución de cada una de ellas y evaluar posibles patrones o comportamientos que puedan influir posteriormente en el comportamiento del modelo.

### Variables Categóricas

```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V1)) +
  geom_bar(fill = "skyblue", color = "black") +
  labs(title = "Distribución de V1", x = "V1", y = "Frecuencia") +
  theme_minimal()
```

-   **Observación:** Podemos observar dos niveles (`a` y `b`), donde **`b` tiene una frecuencia bastante mayor que `a`**. Esto apunta a un **desbalance de clases** que **podría hacer que el modelo favorezca a la categoría mayoritaria**. Este sesgo puede resultar en que **no se generalice correctamente** sobre los casos en `a`, lo que **afectaría a la precisión y la fiabilidad**.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V4)) +
  geom_bar(fill = "skyblue", color = "black") +
  labs(title = "Distribución de V4", x = "V4", y = "Frecuencia") +
  theme_minimal()
```

-   **Observación:** De las tres categorías (`l`, `u`, `y`), **`u` tiene la mayor frecuencia con diferencia** y `l` casi no está representada. La concentración en **`u` sugiere que esta es la categoría representativa**, mientras que las categorías `l` e `y` son casos menos frecuentes o con un comportamiento atípico. **Si `l` es irrelevante o está sesgado, es posible que introduzca ruido** o **imprecisiones en el modelo**, mientras que **`u` podría estar fuertemente correlacionada con la variable objetivo (`V16`)**, opacando a `l` e `y` en la predicción. **Este desbalance podría afectar la precisión del modelo**, pues podría generar **sesgo hacia `u`**.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V5)) +
  geom_bar(fill = "skyblue", color = "black") +
  labs(title = "Distribución de V5", x = "V5", y = "Frecuencia") +
  theme_minimal()
```

-   **Observación:** Distribución parecida a V4 con **predominancia de `g`** y casi nula presencia de `gg`.

```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V6)) +
  geom_bar(fill = "skyblue", color = "black") +
  labs(title = "Distribución de V6", x = "V6", y = "Frecuencia") +
  theme_minimal()
```

-   **Observación:** Presenta múltiples niveles con distribución desigual, donde **`c` es la más dominante con diferencia**, mientras que categorías como **`j` y `r` tienen frecuencias bajas**. La categoría dominante (**`c`**) podría representar **comportamientos más comunes o riesgos identificables**, mientras que las categorías poco comunes (como **`j` y `r`**) podrían ser **casos específicos, erróneos o inusuales**. **La irregularidad de esta variable podría producir ruido a la hora de entrenar el modelo**.
    

```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V7)) +
  geom_bar(fill = "skyblue", color = "black") +
  labs(title = "Distribución de V7", x = "V7", y = "Frecuencia") +
  theme_minimal()
```

-   **Observación:** Similar a **V6**, esta vez podemos ver que la categoría dominante es **`v`**, mientras que variables como **`o`, `z`, `n`, `j` y `dd`** son **minoritarios y casi no tienen frecuencia**. Debido a la diferencia tan clara de frecuencia de las categorías podría generar **imprecisiones a la hora de entrenar el modelo**.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V9)) +
  geom_bar(fill = "skyblue", color = "black") +
  labs(title = "Distribución de V9", x = "V9", y = "Frecuencia") +
  theme_minimal()
```

-   **Observación:** Aquí vemos que la **V9** está compuesta por dos categorías (**`f`** y **`t`**) donde hay prácticamente **igualdad de elementos en cada categoría**. Esta distribución **equitativa** le convierte en candidato a ser una **variable clave a la hora de entrenar el modelo** porque no generaría mucho sesgo hacia ninguna categoría.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V10)) +
  geom_bar(fill = "skyblue", color = "black") +
  labs(title = "Distribución de V10", x = "V10", y = "Frecuencia") +
  theme_minimal()
```

-   **Observación:** Similar a **`V9`**, con niveles **`f`** y **`t`** siendo el primero el **mayoritario**. Esta variable también sigue una **distribución bastante equitativa**, aunque no tanto como **V9**. Podría generar **cierto sesgo**.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V12)) +
  geom_bar(fill = "skyblue", color = "black") +
  labs(title = "Distribución de V12", x = "V12", y = "Frecuencia") +
  theme_minimal()
```

-   **Observación:** Dos niveles (**`f`** y **`t`**), con una **distribución muy similar** y **`f` dominante**. Parecida a las variables **V9** y **V10**.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V13)) +
  geom_bar(fill = "skyblue", color = "black") +
  labs(title = "Distribución de V13", x = "V13", y = "Frecuencia") +
  theme_minimal()
```

-   **Observación:** Tres niveles (`g`, `p`, y `s`), con **`g` como dominante**. Este podría estar asociada con el **comportamiento general**, mientras que las categorías **`p` y `s` presentan una frecuencia minoritaria**, lo que puede hacerlas **menos útiles**, o que sean **valores atípicos**.


### Variables Numéricas

```{r, fig.width=5, fig.height=3}
# Graficar cada variable numérica individualmente
ggplot(credit, aes(x = V2)) +
  geom_histogram(fill = "orange", color = "black", bins = 20, alpha = 0.7) +
  labs(title = "Distribución de V2", x = "V2", y = "Frecuencia") +
  geom_vline(xintercept = mean(credit$V2), col = "blue") +
  geom_vline(xintercept = median(credit$V2), col = "red") +
  geom_rug(data=credit, aes(x=V2, y=0), sides="b", position="jitter") +
  theme_minimal()
```

-   **Observación:** Presenta **asimetría positiva**, con **concentración en valores bajos** y **cola extendida hacia la derecha**. Esto refleja que la **mayoría de los casos son pequeños**, mientras que los **valores altos pueden ser atípicos**.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V3)) +
  geom_histogram(fill = "orange", color = "black", bins = 20, alpha = 0.7) +
  labs(title = "Distribución de V3", x = "V3", y = "Frecuencia") +
  geom_vline(xintercept = mean(credit$V3), col = "blue") +
  geom_vline(xintercept = median(credit$V3), col = "red") +
  geom_rug(data=credit, aes(x=V3, y=0), sides="b", position="jitter") +
  theme_minimal()
```

-   **Observación:** También presenta **asimetría positiva**, aunque con **mayor concentración en valores bajos**. Al igual que con `V2`, los **valores altos pueden representar casos atípicos**.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V8)) +
  geom_histogram(fill = "orange", color = "black", bins = 20, alpha = 0.7) +
  labs(title = "Distribución de V8", x = "V8", y = "Frecuencia") +
  geom_vline(xintercept = mean(credit$V8), col = "blue") +
  geom_vline(xintercept = median(credit$V8), col = "red") +
  geom_rug(data=credit, aes(x=V8, y=0), sides="b", position="jitter") +
  theme_minimal()
```

-   **Observación:** Similar a **`V3`**, lo que sugiere un **comportamiento similar**.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V11)) +
  geom_histogram(fill = "orange", color = "black", bins = 20, alpha = 0.7) +
  labs(title = "Distribución de V11", x = "V11", y = "Frecuencia") +
  geom_vline(xintercept = mean(credit$V11), col = "blue") +
  geom_vline(xintercept = median(credit$V11), col = "red") +
  geom_rug(data=credit, aes(x=V11, y=0), sides="b", position="jitter") +
  theme_minimal()
```

-   **Observación:** Valores **bajos predominantes**, existe una **cola extendida** hacia valores **altos**. Similar a los anteriores.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V14)) +
  geom_histogram(fill = "orange", color = "black", bins = 20, alpha = 0.7) +
  labs(title = "Distribución de V14", x = "V14", y = "Frecuencia") +
  geom_vline(xintercept = mean(credit$V14), col = "blue") +
  geom_vline(xintercept = median(credit$V14), col = "red") +
  geom_rug(data=credit, aes(x=V14, y=0), sides="b", position="jitter") +
  theme_minimal()
```

-   **Observación:** Valores mayoritariamente **bajos**. Similar a **V2**, pero con una **concentración muy alta** al principio.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x = V15)) +
  geom_histogram(fill = "orange", color = "black", bins = 20, alpha = 0.7) +
  labs(title = "Distribución de V15", x = "V15", y = "Frecuencia") +
  geom_vline(xintercept = mean(credit$V15), col = "blue") +
  geom_vline(xintercept = median(credit$V15), col = "red") +
  geom_rug(data=credit, aes(x=V15, y=0), sides="b", position="jitter") +
  theme_minimal()
```

-   **Observación:** Similar a **V3** o **V8**, pero con una **concentración en el rango bajo** mucho más **exagerada** y una **cola extendida muy larga** aunque con **muy poca frecuencia**. Es posible que se traten de **datos atípicos** aunque podría estar relacionado con la **aprobación del crédito**.


## Análisis de cada variable en función del crédito

### Variables Categóricas

```{r, fig.width=5, fig.height=3}
# Análisis de cada variable categórica frente a la aprobación del crédito (V16)
ggplot(data = credit, aes(x = V1, fill = V16)) +
  geom_bar() +
  labs(y = "Proporción", x = 'V1') + ggtitle('V1 vs V16') +
  theme_minimal()
```

Hay una **mayor proporción** de **aprobaciones** en `b`, que en `a`. Se ve una **distribución desigual**, por lo que podríamos pensar que **es más probable que `b` tenga el crédito aprobado**, comparado con `a`.


```{r, fig.width=5, fig.height=3}
ggplot(data = credit, aes(x = V4, fill = V16)) +
  geom_bar() +
  labs(y = "Proporción", x = 'V4') + ggtitle('V4 vs V16') +
  theme_minimal()
```

En este caso, `u` es la que parece estar relacionada con un **mayor porcentaje** de ser **aprobado**, pero a su vez, también con el **mayor porcentaje** de **no** ser **aprobado**. La **categoría `y`** tiene una **predominancia a no ser aprobada**, lo que podría ayudar a predecir.

Por otro lado, `l` tiene muy **poca cantidad de datos**. Podríamos pensar que es **ruido** o un caso muy especial, en el que se aprueba el crédito. **Más adelante**, lo **analizaremos en mayor profundidad** y veremos si los quitamos o no.


```{r, fig.width=5, fig.height=3}
ggplot(data = credit, aes(x = V5, fill = V16)) +
  geom_bar() +
  labs(y = "Proporción", x = 'V5') + ggtitle('V5 vs V16') +
  theme_minimal()
```

En el caso `g`, podemos apreciar una **distribución más homogénea** entre **aprobados** y **no aprobados**. Las demás, son **proporciones más bajas**.

Y al igual que antes, en `gg`, podríamos considerar que son **datos anómalos** o **ruido**. **Debemos analizarlo posteriormente** para ver si los seguimos considerando, o no.


```{r, fig.width=5, fig.height=3}
ggplot(data = credit, aes(x = V6, fill = V16)) +
  geom_bar() +
  labs(y = "Proporción", x = 'V6') + ggtitle('V6 vs V16') +
  theme_minimal()
```

Se ve como la **categoría `c`** es la que **mayor número de no aprobado tiene**, por lo que podría ser de **alto riesgo**, aunque no es el que más tiene en término de porcentaje. Por ejemplo, **`ff`** e **`i`** tienen una **clara tendencia a no ser aprobada**. Por otro lado, el nivel **`cc`** y **`x`** son los que **mayor tasa de aprobado tienen**.

Además, volvemos a encontrarnos con categorías que tienen **pocos datos**, **como `j` y `r`**, en comparación con las demás.

Podríamos pensar que hay **subgrupos de niveles** que podrían ser **relevantes** para nuestra **función objetivo**, pues en general, **V6 muestra diferencias** en la tasa de aprobado **según la categoría**.


```{r, fig.width=5, fig.height=3}
ggplot(data = credit, aes(x = V7, fill = V16)) +
  geom_bar() +
  labs(y = "Proporción", x = 'V7') + ggtitle('V7 vs V16') +
  theme_minimal()
```

Vemos que **`ff`** y **`v`** son de las **categorías** con **mayor probabilidad de no tener el crédito aprobado**, **`h`** es una **categoría bastante favorable a tener el crédito aprobado** y tenemos categorías con **pocos datos**, como `dd`, `j`, `n`, `o` y `z`, que podrían ser **valores atípicos**.

Es así, que **V7 también parece ser una variable que afecta a la probabilidad** de tener o no el aprobado.


```{r, fig.width=5, fig.height=3}
ggplot(data = credit, aes(x = V9, fill = V16)) +
  geom_bar() +
  labs(y = "Proporción", x = 'V9') + ggtitle('V9 vs V16') +
  theme_minimal()
```

La categoría **`f`** es **muy favorable a que no le aprueben el crédito**, mientras que **`t` sí tiene una mayor tasa de aprobado**. Esto, sumado a que tienen una **frecuencia muy parecida**, hace que sea una variable **muy a tener en cuenta a la hora de clasificar**.

En conclusión, **V9 es una variable bastante discriminante**.


```{r, fig.width=5, fig.height=3}
ggplot(data = credit, aes(x = V10, fill = V16)) +
  geom_bar() +
  labs(y = "Proporción", x = 'V10') + ggtitle('V10 vs V16') +
  theme_minimal()
```

En este caso, **`f`** es **más favorable a que no le aprueben el crédito**, mientras que **`t`** **se relaciona con un porcentaje mayor de aprobados**.

**V10 podría ser una variable relevante para predecir el aprobado.**


```{r, fig.width=5, fig.height=3}
ggplot(data = credit, aes(x = V12, fill = V16)) +
  geom_bar() +
  labs(y = "Proporción", x = 'V12') + ggtitle('V12 vs V16') +
  theme_minimal()
```

La diferencia entre ambas categorías **no parece ser muy grande**. De modo que **parece que no hay relación entre ellas**.


```{r, fig.width=5, fig.height=3}
ggplot(data = credit, aes(x = V13, fill = V16)) +
  geom_bar() +
  labs(y = "Proporción", x = 'V13') + ggtitle('V13 vs V16') +
  theme_minimal()
```

La categoría **`g`** es **bastante dominante**, donde **predominan los rechazos aunque con poca fuerza**. Mientras que **`p`** y **`s`** tienen **muy pocos datos**. Podrían tratarse de **valores atípicos**.

Esto hace que **V13 no sea potencialmente una buena variable para discriminar**.


### Variables Numéricas

```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x= V16, y= V2, fill= V16)) +
  geom_boxplot() +
  labs(y = "V2", x = 'V16') + ggtitle('V2 vs V16') +
  scale_fill_brewer(palette = "Set2")
```

Se aprecia como la **mediana** de los que tuvieron el crédito **rechazado** es **mayor** que la mediana de los que fueron **aprobados**. Pero no parece haber una **correlación fuerte**.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x= V16, y= V3, fill= V16)) +
  geom_boxplot() +
  labs(y = "V3", x = 'V16') + ggtitle('V3 vs V16') +
  scale_fill_brewer(palette = "Set2")
```

La **mediana** de los **aprobados** es **mayor** que la de los **rechazados**. Por tanto, los valores más altos de **V3** podrían tener mayor posibilidad de ser **aprobados**. Los **rechazados** tienden a agruparse en los **valores bajos**. Parece que **V3** tiene cierta capacidad **discriminante**.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x= V16, y= V8, fill= V16)) +
  geom_boxplot() +
  labs(y = "V8", x = 'V16') + ggtitle('V8 vs V16') +
  scale_fill_brewer(palette = "Set2")
```

Al igual que **V3**, se ve como la **mediana** de los **aprobados** es **mayor** que la de los **rechazados**. Por lo que para **valores más altos**, hay **mayor tasa de aprobados**. Y que los **no aprobados**, se concentran en **valores bajos**. Aquí sí que se ve una **gran diferencia** entre ambos, por lo que quizás sería una buena consideración para el modelo.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x= V16, y= V11, fill= V16)) +
  geom_boxplot() +
  labs(y = "V11", x = 'V16') + ggtitle('V11 vs V16') +
  scale_fill_brewer(palette = "Set2")
```

Con **V11**, volvemos a ver una **relación positiva**, y bastante **discriminante**. Para los **valores mayores, se relacionan principalmente con aprobados**, mientras que los **valores más bajos**, se relacionan con el **rechazo**.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x= V16, y= V14, fill= V16)) +
  geom_boxplot() +
  labs(y = "V14", x = 'V16') + ggtitle('V14 vs V16') +
  scale_fill_brewer(palette = "Set2")
```

La **mediana** de los **no aprobados** es **menor**, concentrándose en los **valores más bajos**, y la **mediana** de los **aprobados** es **mayor**, por lo que la **tasa de aprobados** se asocia a **valores mayores** de **V14**.

Existen **valores atípicos**, pues hay casos que se van hasta **2000**. Más adelante nos encargaremos de **analizarlos**.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x= V16, y= V15, fill= V16)) +
  geom_boxplot() +
  labs(y = "V15", x = 'V16') + ggtitle('V15 vs V16') +
  scale_fill_brewer(palette = "Set2")
```

En este último **boxplot** comparando con **V15**, se ve que tenemos **valores atípicos**, que nos hacen ver la gráfica con **poca claridad**. Más adelante, nos encargaremos de **analizar los valores atípicos**, pero por ahora, solo vamos a **escalar la gráfica** con tal de verlo mejor.


```{r, fig.width=5, fig.height=3}
ggplot(credit, aes(x= V16, y= V15, fill= V16)) +
  geom_boxplot() +
  labs(y = "V15", x = 'V16') + ggtitle('V15 vs V16') +
  scale_fill_brewer(palette = "Set2") +
  coord_cartesian(ylim=c(0, 2000))
```

Podemos ver una **fuerte relación** entre **V15** y la **variable objetivo**. Parece ser que es una **variable altamente discriminante**, y que existe una **correlación positiva**.

## Análisis de Outliers

Con las visualizaciones anteriores, hemos visto que hay datos
considerados atípicos, y que nos encargaríamos de analizarlos. Es
importante tratarlos de cara a los modelos que vayamos a entrenar, pues
**podría influir de mala manera en los resultados**. Es por ello, que los trataremos
con el conjunto de train.

Vamos a pasar las variables numéricas a formato largo, con tal de poder
generar los boxplots para todas, y así visualizarlas de manera conjunta
y detectar valores atípicos.

```{r}
# Transformar el dataset al formato largo
credit_melt <- melt(credit.Datos.Train, measure.vars = which(numeric_columns))
head(credit_melt)

```

```{r, fig.width=5, fig.height=3}

# Boxplots de todas las variables numéricas
ggplot(credit_melt, aes(x = variable, y = value)) +
  geom_boxplot(fill = "skyblue", outlier.color = "red", outlier.shape = 16) +
  labs(
    title = "Boxplots de Variables Numéricas",
    x = "Variables",
    y = "Valores"
  ) +
  theme_minimal()


```

Como podemos ver, **V15** es la **variable** que **más se dispara con diferencia**.

### Boxplots desgregados por Clase (V16)

Si ahora los visualizamos desgregados por la variable objetivo,
podríamos ver cómo los outliers podrían afectar a esta.

```{r, fig.width=5, fig.height=3}
ggplot(credit_melt, aes(x = variable, y = value, fill = V16)) +
  geom_boxplot(outlier.color = "red", outlier.shape = 16) +
  labs(
    title = "Boxplots de Variables Numéricas por Clase",
    x = "Variables",
    y = "Valores",
    fill = "Clase"
  ) +
  theme_minimal()


```

Aquí tambien hay un **número significativo de outliers en V15**,concentrándose en **valores extremos muy altos**.

### Criterio de Turkey

La manera más fácil de tratar outliers es con el **Criterio de Turkey**. Consiste en eliminar los **valores atípicos**, utilizando el **rango intercuartílico (IQR)** para detectar los **límites superior e inferior**. Los que queden fuera de esos límites, se consideran **outliers**.


```{r}
# Función para detectar outliers
detect_outliers <- function(x) {
  Q1 <- quantile(x, 0.25, na.rm = TRUE)
  Q3 <- quantile(x, 0.75, na.rm = TRUE)
  IQR <- Q3 - Q1
  lower_bound <- Q1 - 1.5 * IQR
  upper_bound <- Q3 + 1.5 * IQR
  return(x < lower_bound | x > upper_bound)
}

# Detectar outliers por variable
outliers <- lapply(credit.Datos.Train[numeric_columns], detect_outliers)

# Contar outliers por variable
sapply(outliers, sum)

total_outliers <- sum(sapply(outliers, sum))
total_outliers
```

Aparecen demasiados valores considerados outliers.

### Distancias de Mahalanobis

Vamos a probar con otro criterio, pero esta vez **tiene en cuenta outliers multivariados**. O sea, se consideran las **relaciones entre las distintas variables** y nos ayuda a detectar **patrones anómalos** al patrón general de los datos.


```{r, fig.width=5, fig.height=3}
# Cálculo y visualización de distancias de Mahalanobis

#Calcular vector de medias y matriz de covarianza
mean_vector <- colMeans(credit.Datos.Train[, numeric_columns], na.rm = TRUE)
cov_matrix <- cov(credit.Datos.Train[, numeric_columns], use = "complete.obs")

# Calcular distancias de Mahalanobis
mahalanobis_distances <- mahalanobis(credit.Datos.Train[, numeric_columns], center = mean_vector, cov = cov_matrix)

# Umbral para identificar outliers (97.5% de confianza)
threshold <- qchisq(0.975, df = ncol(credit.Datos.Train[, numeric_columns]))

# Identificar outliers
outliers <- which(mahalanobis_distances > threshold)

plot(mahalanobis_distances,
     main = "Distancias de Mahalanobis",
     xlab = "Índice de Observación",
     ylab = "Distancia de Mahalanobis",
     pch = 19, col = "blue")
abline(h = threshold, col = "red", lty = 2) # Línea roja para el umbral

# outliers totales
num_outliers <- length(outliers)
cat("Número de outliers identificados:", num_outliers, "\n")

```

Comparado con el criterio de Turkey, hay muchos menos outliers.

### Acciones a tomar ante los outliers

Como no tenemos ningún contexto de lo que es cada variable, pues están codificadas, y no tenemos a ningún experto que nos pudiera ayudar, vamos a optar por **quitar los outliers**. Además, vamos a emplear el **criterio de las distancias de Mahalanobis**, pues con el **criterio de Turkey**, se elimina una **gran cantidad de datos**, y tampoco nos podemos permitir eso, al ser un dataset bastante **"pequeño"** ya de por sí. Por no decir que es **más robusto** un **análisis multivariado** que tiene en cuenta las **relaciones entre las variables**.


```{r}

# Verificar el tamaño del dataset antes de eliminar los outliers
cat("Número de filas antes de eliminar outliers:", nrow(credit.Datos.Train), "\n")

# Eliminar los outliers del dataset
credit.Datos.Train <- credit.Datos.Train[-outliers, ]

# Verificar el tamaño del dataset después de eliminar los outliers
cat("Número de filas después de eliminar outliers:", nrow(credit.Datos.Train), "\n")
```

### Revisualizar los datos después de quitar los outliers

```{r graficos-grid2, fig.width=10, fig.height=12, echo=FALSE, warning=FALSE}

# Crear lista de gráficos para variables numéricas
numeric_plots <- lapply(numeric_columns2, function(col) {
  ggplot(credit.Datos.Train, aes_string(y = col)) +
    geom_boxplot(fill = "orange", color = "black", outlier.color = "red") +
    labs(title = paste("Boxplot de", col), y = col, x = "") +
    theme_minimal()
})

# Mostrar gráficos numéricos en un grid
grid.arrange(grobs = numeric_plots, ncol = 2, top = "Variables Numéricas")


```

Los **boxplots** han cambiado con respecto a los que vimos anteriormente. Pero se observa que aún así, en **V15 sigue teniendo una gran cantidad de valores extremadamente elevados**. Como hemos empleado un **enfoque multivariado** para eliminar los **outliers**, esto sugiere que esos **valores tan extremos no son inconsistentes con el resto de variables del dataset**. Esto podría indicar que **V15** refleja algún tipo de **métrica económica**, o algo relacionado con el **dinero** de las personas, lo que justificaría la presencia de valores altos en algunos casos.

No obstante, más adelante, nos encargaremos de **normalizar los datos**, con tal de que esta diferencia de escala no afecte a los modelos. Si tuviéramos **contexto** de esta variable, sería más fácil decidir si recortar de una manera más agresiva o dejarla así, por tener un significado útil. Pero **nos fiaremos del enfoque multivariado**.



## Comprobar distribución

Es interesante comprobar las **distribuciones** que siguen las variables, **a la hora de aplicar transformaciones adicionales**, pues **hay modelos que necesitan tener los datos normalizados** o que sigan una cierta forma. Para ello, vamos a realizar los **Q-Q plots** correspondientes a las distribuciones más comunes, como la **normal**, **binomial**, **exponencial**, **geométrica** o **poisson**.


### Distribución Normal

Las variables numericas, a simple vista con las gráficas anteriores, 
no parecen seguir una distribución normal. Sino, una geometrica o exponencial.

De todos modos, vamos a mirar si las variables se adaptan a una
distribución normal.

```{r graficas Q-Q, fig.width=14, fig.height=10}

p1 = ggplot(data=credit.Datos.Train, aes(sample=credit.Datos.Train[,2])) +
  ggtitle("QQ plot V2") +
  geom_qq() + 
  stat_qq_line() + 
  xlab("Distribución teórica") + ylab("Distribución muestral")

p2 = ggplot(data=credit.Datos.Train, aes(sample=credit.Datos.Train[,3])) +
  ggtitle("QQ plot V3") +
  geom_qq() + 
  stat_qq_line() + 
  xlab("Distribución teórica") + ylab("Distribución muestral")

p3 = ggplot(data=credit.Datos.Train, aes(sample=credit.Datos.Train[,8])) +
  ggtitle("QQ plot V8") +
  geom_qq() + 
  stat_qq_line() + 
  xlab("Distribución teórica") + ylab("Distribución muestral")

p4 = ggplot(data=credit.Datos.Train, aes(sample=credit.Datos.Train[,11])) +
  ggtitle("QQ plot V11") +
  geom_qq() + 
  stat_qq_line() + 
  xlab("Distribución teórica") + ylab("Distribución muestral")

p5 = ggplot(data=credit.Datos.Train, aes(sample=credit.Datos.Train[,14])) +
  ggtitle("QQ plot V14") +
  geom_qq() + 
  stat_qq_line() + 
  xlab("Distribución teórica") + ylab("Distribución muestral")

p6 = ggplot(data=credit.Datos.Train, aes(sample=credit.Datos.Train[,15])) +
  ggtitle("QQ plot V15") +
  geom_qq() + 
  stat_qq_line() + 
  xlab("Distribución teórica") + ylab("Distribución muestral")

#Organizando los gráficos en un solo lienzo
grid.arrange(p1, p2, p3, p4, p5, p6, nrow=2)


```

En las QQ plots ajustadas a la distribución normal, se puede observar que
**ninguna de las variables sigue un comportamiento cercano a la normalidad**. 
Además, las variables muestran colas alargadas (asimetría positiva) en 
los valores altos.

### Distribución Binomial

```{r graficas Q-Q Binomial, fig.width=14, fig.height=10}

# Establecer los parámetros para cada variable
# Usamos como número de ensayos (n) la máxima observación + 1, por simplicidad
n_v2 <- max(credit.Datos.Train[, 2]) + 1
n_v3 <- max(credit.Datos.Train[, 3]) + 1
n_v8 <- max(credit.Datos.Train[, 8]) + 1
n_v11 <- max(credit.Datos.Train[, 11]) + 1
n_v14 <- max(credit.Datos.Train[, 14]) + 1
n_v15 <- max(credit.Datos.Train[, 15]) + 1

# Probabilidad de éxito (p) como media / n
p_binom_v2 <- mean(credit.Datos.Train[, 2]) / n_v2
p_binom_v3 <- mean(credit.Datos.Train[, 3]) / n_v3
p_binom_v8 <- mean(credit.Datos.Train[, 8]) / n_v8
p_binom_v11 <- mean(credit.Datos.Train[, 11]) / n_v11
p_binom_v14 <- mean(credit.Datos.Train[, 14]) / n_v14
p_binom_v15 <- mean(credit.Datos.Train[, 15]) / n_v15

# Crear gráficos Q-Q para cada variable
p1_binomial <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 2])) +
  stat_qq(distribution = stats::qbinom, dparams = list(size = n_v2, prob = p_binom_v2)) +
  stat_qq_line(distribution = stats::qbinom, dparams = list(size = n_v2, prob = p_binom_v2)) +
  ggtitle("QQ plot - Distribución Binomial (V2)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p2_binomial <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 3])) +
  stat_qq(distribution = stats::qbinom, dparams = list(size = n_v3, prob = p_binom_v3)) +
  stat_qq_line(distribution = stats::qbinom, dparams = list(size = n_v3, prob = p_binom_v3)) +
  ggtitle("QQ plot - Distribución Binomial (V3)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p3_binomial <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 8])) +
  stat_qq(distribution = stats::qbinom, dparams = list(size = n_v8, prob = p_binom_v8)) +
  stat_qq_line(distribution = stats::qbinom, dparams = list(size = n_v8, prob = p_binom_v8)) +
  ggtitle("QQ plot - Distribución Binomial (V8)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p4_binomial <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 11])) +
  stat_qq(distribution = stats::qbinom, dparams = list(size = n_v11, prob = p_binom_v11)) +
  stat_qq_line(distribution = stats::qbinom, dparams = list(size = n_v11, prob = p_binom_v11)) +
  ggtitle("QQ plot - Distribución Binomial (V11)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p5_binomial <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 14])) +
  stat_qq(distribution = stats::qbinom, dparams = list(size = n_v14, prob = p_binom_v14)) +
  stat_qq_line(distribution = stats::qbinom, dparams = list(size = n_v14, prob = p_binom_v14)) +
  ggtitle("QQ plot - Distribución Binomial (V14)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p6_binomial <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 15])) +
  stat_qq(distribution = stats::qbinom, dparams = list(size = n_v15, prob = p_binom_v15)) +
  stat_qq_line(distribution = stats::qbinom, dparams = list(size = n_v15, prob = p_binom_v15)) +
  ggtitle("QQ plot - Distribución Binomial (V15)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

# Organizar gráficos en un grid
grid.arrange(p1_binomial, p2_binomial, p3_binomial, p4_binomial, p5_binomial, p6_binomial, nrow = 2)


```

Observando las gráficas, podemos ver que **tampoco se ajusta a una distribución binomial**. Especialmente en variables como V3 y V11, donde los valores esperados no coinciden con los observados.

### Distribución Exponencial

```{r graficas Q-Q Exponencial, fig.width=14, fig.height=10}

# Calcular el parámetro lambda (rate) para cada variable
lambda_exp_v2 <- 1 / mean(credit.Datos.Train[, 2])
lambda_exp_v3 <- 1 / mean(credit.Datos.Train[, 3])
lambda_exp_v8 <- 1 / mean(credit.Datos.Train[, 8])
lambda_exp_v11 <- 1 / mean(credit.Datos.Train[, 11])
lambda_exp_v14 <- 1 / mean(credit.Datos.Train[, 14])
lambda_exp_v15 <- 1 / mean(credit.Datos.Train[, 15])

# Crear gráficos Q-Q para cada variable
p1_exponential <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 2])) +
  stat_qq(distribution = stats::qexp, dparams = list(rate = lambda_exp_v2)) +
  stat_qq_line(distribution = stats::qexp, dparams = list(rate = lambda_exp_v2)) +
  ggtitle("QQ plot - Distribución Exponencial (V2)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p2_exponential <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 3])) +
  stat_qq(distribution = stats::qexp, dparams = list(rate = lambda_exp_v3)) +
  stat_qq_line(distribution = stats::qexp, dparams = list(rate = lambda_exp_v3)) +
  ggtitle("QQ plot - Distribución Exponencial (V3)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p3_exponential <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 8])) +
  stat_qq(distribution = stats::qexp, dparams = list(rate = lambda_exp_v8)) +
  stat_qq_line(distribution = stats::qexp, dparams = list(rate = lambda_exp_v8)) +
  ggtitle("QQ plot - Distribución Exponencial (V8)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p4_exponential <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 11])) +
  stat_qq(distribution = stats::qexp, dparams = list(rate = lambda_exp_v11)) +
  stat_qq_line(distribution = stats::qexp, dparams = list(rate = lambda_exp_v11)) +
  ggtitle("QQ plot - Distribución Exponencial (V11)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p5_exponential <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 14])) +
  stat_qq(distribution = stats::qexp, dparams = list(rate = lambda_exp_v14)) +
  stat_qq_line(distribution = stats::qexp, dparams = list(rate = lambda_exp_v14)) +
  ggtitle("QQ plot - Distribución Exponencial (V14)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p6_exponential <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 15])) +
  stat_qq(distribution = stats::qexp, dparams = list(rate = lambda_exp_v15)) +
  stat_qq_line(distribution = stats::qexp, dparams = list(rate = lambda_exp_v15)) +
  ggtitle("QQ plot - Distribución Exponencial (V15)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

# Organizar gráficos en un grid
grid.arrange(p1_exponential, p2_exponential, p3_exponential, 
             p4_exponential, p5_exponential, p6_exponential, nrow = 2)

```

Viedno estas gráficas, la **única que se ajusta razonadamente a la distribución exponencial es la V8**, especialmente en los valores más bajos, mientras que 
las demás variables, como **V2** y **V15**, **presentan grandes desviaciones en las colas**. 

Esto indica que V8 podría modelar procesos decayentes, como tiempos.

### Distribución Geométrica

```{r graficas Q-Q Geométrica, fig.width=14, fig.height=10}

# Calcular el parámetro p para cada variable
p_geom_v2 <- 1 / mean(credit.Datos.Train[, 2])
p_geom_v3 <- 1 / mean(credit.Datos.Train[, 3])
p_geom_v8 <- 1 / mean(credit.Datos.Train[, 8])
p_geom_v11 <- 1 / mean(credit.Datos.Train[, 11])
p_geom_v14 <- 1 / mean(credit.Datos.Train[, 14])
p_geom_v15 <- 1 / mean(credit.Datos.Train[, 15])

# Crear gráficos Q-Q para cada variable
p1_geometric <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 2])) +
  stat_qq(distribution = stats::qgeom, dparams = list(prob = p_geom_v2)) +
  stat_qq_line(distribution = stats::qgeom, dparams = list(prob = p_geom_v2)) +
  ggtitle("QQ plot - Distribución Geométrica (V2)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p2_geometric <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 3])) +
  stat_qq(distribution = stats::qgeom, dparams = list(prob = p_geom_v3)) +
  stat_qq_line(distribution = stats::qgeom, dparams = list(prob = p_geom_v3)) +
  ggtitle("QQ plot - Distribución Geométrica (V3)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p3_geometric <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 8])) +
  stat_qq(distribution = stats::qgeom, dparams = list(prob = p_geom_v8)) +
  stat_qq_line(distribution = stats::qgeom, dparams = list(prob = p_geom_v8)) +
  ggtitle("QQ plot - Distribución Geométrica (V8)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p4_geometric <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 11])) +
  stat_qq(distribution = stats::qgeom, dparams = list(prob = p_geom_v11)) +
  stat_qq_line(distribution = stats::qgeom, dparams = list(prob = p_geom_v11)) +
  ggtitle("QQ plot - Distribución Geométrica (V11)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p5_geometric <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 14])) +
  stat_qq(distribution = stats::qgeom, dparams = list(prob = p_geom_v14)) +
  stat_qq_line(distribution = stats::qgeom, dparams = list(prob = p_geom_v14)) +
  ggtitle("QQ plot - Distribución Geométrica (V14)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p6_geometric <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 15])) +
  stat_qq(distribution = stats::qgeom, dparams = list(prob = p_geom_v15)) +
  stat_qq_line(distribution = stats::qgeom, dparams = list(prob = p_geom_v15)) +
  ggtitle("QQ plot - Distribución Geométrica (V15)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

# Organizar gráficos en un grid
grid.arrange(p1_geometric, p2_geometric, p3_geometric, 
             p4_geometric, p5_geometric, p6_geometric, nrow = 2)

```

Una vez más no parece que ninguna se ajuste a una distribución geométrica tampoco.

### Distribución Poisson

```{r graficas Q-Q Poisson, fig.width=14, fig.height=10}

# Calcular el parámetro lambda (la media de cada variable)
lambda_pois_v2 <- mean(credit.Datos.Train[, 2])
lambda_pois_v3 <- mean(credit.Datos.Train[, 3])
lambda_pois_v8 <- mean(credit.Datos.Train[, 8])
lambda_pois_v11 <- mean(credit.Datos.Train[, 11])
lambda_pois_v14 <- mean(credit.Datos.Train[, 14])
lambda_pois_v15 <- mean(credit.Datos.Train[, 15])

# Crear gráficos Q-Q para cada variable
p1_poisson <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 2])) +
  stat_qq(distribution = stats::qpois, dparams = list(lambda = lambda_pois_v2)) +
  stat_qq_line(distribution = stats::qpois, dparams = list(lambda = lambda_pois_v2)) +
  ggtitle("QQ plot - Distribución Poisson (V2)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p2_poisson <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 3])) +
  stat_qq(distribution = stats::qpois, dparams = list(lambda = lambda_pois_v3)) +
  stat_qq_line(distribution = stats::qpois, dparams = list(lambda = lambda_pois_v3)) +
  ggtitle("QQ plot - Distribución Poisson (V3)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p3_poisson <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 8])) +
  stat_qq(distribution = stats::qpois, dparams = list(lambda = lambda_pois_v8)) +
  stat_qq_line(distribution = stats::qpois, dparams = list(lambda = lambda_pois_v8)) +
  ggtitle("QQ plot - Distribución Poisson (V8)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p4_poisson <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 11])) +
  stat_qq(distribution = stats::qpois, dparams = list(lambda = lambda_pois_v11)) +
  stat_qq_line(distribution = stats::qpois, dparams = list(lambda = lambda_pois_v11)) +
  ggtitle("QQ plot - Distribución Poisson (V11)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p5_poisson <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 14])) +
  stat_qq(distribution = stats::qpois, dparams = list(lambda = lambda_pois_v14)) +
  stat_qq_line(distribution = stats::qpois, dparams = list(lambda = lambda_pois_v14)) +
  ggtitle("QQ plot - Distribución Poisson (V14)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

p6_poisson <- ggplot(data = credit.Datos.Train, aes(sample = credit.Datos.Train[, 15])) +
  stat_qq(distribution = stats::qpois, dparams = list(lambda = lambda_pois_v15)) +
  stat_qq_line(distribution = stats::qpois, dparams = list(lambda = lambda_pois_v15)) +
  ggtitle("QQ plot - Distribución Poisson (V15)") +
  xlab("Cuantiles teóricos") + ylab("Cuantiles observados") +
  theme_minimal()

# Organizar gráficos en un grid
grid.arrange(p1_poisson, p2_poisson, p3_poisson, p4_poisson, p5_poisson, p6_poisson, nrow = 2)

```

Tampoco parece que sigan una distribución de Poisson, pues también se muestran 
grandes desviaciones en las colas, en los valores extremos.

En **conclusión**, **ninguna** de las variables **parece ajustarse completamente a una distribución concreta**, salvo **V8** en la **distribución exponencial**. Es por ello, que **a la hora de entrenar debemos considerar la normalización** de estos, o usar **escalados**. Para el **PCA** y el **modelo de redes neuronales** es especialmente importante este paso.


## Análisis multivariable

El **análisis de correlación** nos permite identificar **relaciones lineales** entre las diferentes **variables numéricas** del dataset. Este paso es importante para entender cómo interactúan las variables entre sí y para detectar **patrones relevantes** que puedan influir a la hora de entrenar los modelos. **Para ello, vamos a usar todo el dataset**, para darnos una idea general de las **relaciones entre las variables**.


### Visualización de la Matriz de Correlación

```{r}
corrplot::corrplot(cor(credit[, numeric_columns]))
```

Aquí podemos ver que generalmente hay **poca correlación**, lo que indica que son **independientes**. Aunque parece que hay una **relación moderada entre las parejas** **V2-V8**, **V8-V11** y **V11-V15**. Por ejemplo, **V2 y V8** están **moderadamente correlacionadas de forma positiva**, lo que sugiere que cuando **V2 aumenta**, **V8** también tiende a hacerlo.

### Visualización Gráfica de Pares

Aunque ya hemos visto las relaciones lineales que hay, vamos a emplear también la **gráfica de pares**, pues nos ofrece un **análisis más completo** en el que podemos
ver las **dispersiones**, y la **inclusión según la variable objetivo**. Nos va a servir para complementar la matriz de correlaciones.

```{r}
ggpairs(
  credit[, numeric_columns],
  aes(color = credit$V16), # variable de salida
  title = "Relación entre atributos numéricos en Credit Dataset"
)
```

- **Distribuciones:** Las variables **V2**, **V3**, **V8**, **V11**, **V14** y **V15** muestran **asimetrías fuertes**, con concentración en **valores bajos**. Por esto, hacer **transformaciones de escalado** o **normalización** es una buena opción.

- **Relación con V16:** Las clases **aprobadas** tienden a concentrarse en **valores altos** de **V8** y **V15**, por lo que estas variables son **bastante discriminantes**.

- **Relaciones no lineales:** Algunas combinaciones, como **V8 y V11**, muestran **patrones no lineales**, útiles para modelos como **Random Forest** o **Gradient Boosting**.

En general, los **valores bajos** tienden a estar asociados con **créditos no aprobados**, mientras que los **valores altos**, aunque menos comunes, parecen ser **determinantes importantes** para la **aprobación del crédito**. Además, como hemos visto anteriormente, la gráfica de pares confirma que las parejas **V2-V8**, **V8-V11** y **V11-V15** muestran una **correlación moderada**, y que en el resto la mayoría de correlaciones son **bajas**.

## Análisis de Componentes Principales (PCA)
El PCA permite reducir la dimensionalidad de los datos combinando las variables 
numéricas y categóricas (transformadas en dummies) en un análisis único. 
De esta manera, capturamos las relaciones conjuntas entre todas las variables.

### Preparación de los Datos
Primero, convertimos las variables categóricas en variables dummy, y las combinamos 
con las numéricas.

```{r}
#Separar la variable objetivo (V16) antes de procesar
#Las transformaciones realizadas al conjunto de entrenamiento
#tambien se deben realizar al conjunto de tests

credit_data_train <- credit.Datos.Train[, -which(names(credit) == "V16")]
credit_data_test <- credit.Datos.Test[, -which(names(credit) == "V16")]

# Convertir variables categóricas a dummies
categorical_columns_no16 <- names(credit_data_train)[sapply(credit_data_train, is.factor)]

#Solo se realiza sobre el conjunto de entrenamiento
dummies_modelo <- dummyVars(~ ., data = credit_data_train[, categorical_columns_no16, 
                  drop = FALSE], fullRank = TRUE)

#Se transforman ambos conjuntos
credit_dummies_train <- predict(dummies_modelo, newdata = credit_data_train[, 
                        categorical_columns_no16, drop = FALSE])
credit_dummies_test <- predict(dummies_modelo, newdata = credit_data_test[, 
                       categorical_columns_no16, drop = FALSE])

# Combinar variables numéricas con las dummies
credit_dummies_train <- cbind(credit_data_train[, numeric_columns2, drop = FALSE], 
                    credit_dummies_train)
credit_dummies_test <- cbind(credit_data_test[, numeric_columns2, drop = FALSE],
                   credit_dummies_test)

```


### Aplicación del PCA
El PCA se aplica directamente al dataset combinado y estandarizado, poniendo
scale = TRUE. 
A continuación, evaluamos la proporción de varianza explicada por 
cada componente y graficamos los resultados.

```{r}
# Aplicar PCA 
pca_result <- prcomp(credit_dummies_train, scale = TRUE)

# Resumen del PCA
summary(pca_result)

# Calcular proporción de varianza explicada y acumulada
explained_variance <- pca_result$sdev^2 / sum(pca_result$sdev^2)
cumulative_variance <- cumsum(explained_variance)

biplot(pca_result, scale = 0,cex=0.7)

```

Viendo el **summary(pca_result)**, podemos decir lo siguiente:

-   **PC1**: Explica el **9.5%** de la varianza, siendo el componente más
    relevante.
-   **PC2**: Explica un **6.88%**, acumulando un **16.39%** de la varianza explicada
    entre los dos.
-   **PC3**: Explica un **6.44%**, llevando el total acumulado al **22.83%**.

Estos primeros componentes principales son los que **capturan la mayor parte de la estructura y patrones** del dataset.

Los **componentes a partir del PC16 en adelante** explican una **proporción muy pequeña de varianza** y cada vez más insignificante (prácticamente **menos del 2%**). Y los **últimos componentes** (**PC35 a PC37**) prácticamente **no contribuyen a la varianza**.

Por otro lado, el **biplot** muestra las **variables originales** y las **observaciones proyectadas** en los primeros dos componentes (**PC1 y PC2**). Las **flechas más largas** representan **variables con mayor impacto** en la varianza, como **V4.u**. Las **agrupaciones de puntos** podrían indicarnos **patrones o similitudes** entre los datos.


```{r, fig.width=5, fig.height=3}
scores <- as.data.frame(pca_result$x)
ggplot(scores, aes(x = PC1, y = PC2)) +
  geom_point() +
  labs(title = "Proyección de puntos",
       x = "PC1", y = "PC2") +
  theme_minimal()

```

Con este gráfico podemos apreciar mejor los puntos, que con el biplot.
Podemos ver que parece haber 2 grupos principales. El grupo de la
izquierda parece estar más agrupado, mientras que el de la derecha
muestra una mayor dispersión.

### Selección de las componentes principales

Determinamos el número de componentes principales necesarios para
retener al menos el **90%** de la **varianza acumulada**.

```{r}
# Determinar el número de componentes para retener el 90% de varianza
num_components <- which(cumulative_variance >= 0.9)[1]
num_components
```

De acuerdo con esto, **los primeros 24 componentes acumulan más del 90% dela varianza total**. Esto significa que el dataset puede ser proyectado a un espacio de 24 dimensiones, en lugar de las originales, **sin perder demasiada información**.

### Visualización de la varianza explicada

```{r, fig.width=10, fig.height=6}

# Gráfico de varianza explicada
explvar_plot <- ggplot(data = data.frame(Component = seq_along(explained_variance),
                                         Variance = explained_variance),
                       aes(x = Component, y = Variance)) +
  geom_line() + geom_point() +
  labs(title = "Varianza Explicada por Componentes Principales",
       x = "Componente Principal", y = "Proporción de Varianza Explicada") +
  theme_minimal()

# Gráfico de varianza acumulada
cumvar_plot <- ggplot(data = data.frame(Component = seq_along(cumulative_variance),
                                        CumulativeVariance = cumulative_variance),
                      aes(x = Component, y = CumulativeVariance)) +
  geom_line() + geom_point() +
  labs(title = "Varianza Acumulada por Componentes Principales",
       x = "Componente Principal", y = "Varianza Acumulada") +
  theme_minimal()

# Mostrar gráficos lado a lado
grid.arrange(explvar_plot, cumvar_plot, ncol = 2)
```

El gráfico de proporción de varianza explicada muestra una curva decreciente, indicando que **los primeros componentes son los más importantes**.

Por otro lado, la curva de varianza acumulada muestra que el **50% de la varianza se alcanza con los primeros 9 componentes**. Y que el **90% de la varianza se acumula con los primeros 24 componentes**.

### Examinar los loadings

**Los loadings nos muestran cómo las variables originales contribuyen a cada componente principal**. Los valores absolutos más altos, indican las variables con mayor influencia. Podemos extraerlos y analizar los más importantes:

```{r}
# Extraer los loadings
loadings <- pca_result$rotation

# Mostrar los loadings
head(loadings)

#variables más influyentes para los primeros componentes
apply(loadings[, 1:num_components], 2, function(x) {
  abs_x <- abs(x)
  sorted_indices <- order(abs_x, decreasing = TRUE) # Ordenar por importancia
  names(abs_x)[sorted_indices[1:5]]            # 5 primeras mas importantes   
})
```

Con esta función se puede ver la importancia, por ejemplo para PC1 la variable más importante es V4.y, y así se puede mirar para las demás componentes.


### Proyección al nuevo espacio
Usamos las componentes principales seleccionadas para proyectar los datos al espacio reducido.

```{r}

# Transformamos los datos con el pca anteriormente calculado
credit.Datos.Train.Transf.PCA <- predict(pca_result, credit_dummies_train)
credit.Datos.Test.Transf.PCA <- predict(pca_result, credit_dummies_test)

#Nos quedamos solo con los componentes principales que describen el 90% de la varianza
credit.Datos.Train.Transf.PCA <- credit.Datos.Train.Transf.PCA[,1:num_components]
credit.Datos.Test.Transf.PCA <- credit.Datos.Test.Transf.PCA[,1:num_components]

# Ver los resultados
head(credit.Datos.Train.Transf.PCA)
head(credit.Datos.Test.Transf.PCA)


```

Básicamente, el **PCA nos ha permitido reducir la dimensionalidad del dataset**, **reteniendo el 90% de la varianza acumulada** con solo **24 componentes principales**. Así, podríamos mejorar la eficiencia y eliminar redundancias en los datos, a la hora de entrenar modelos.

# Modelos Predictivos

## Normalización de los datos

En los boxplots anteriores, hemos visto como hay escalas diferentes. Por
ejemplo, V15 va desde 0 hasta 500, y V14 va de 100 a alrededor de 250.
Por ello, vamos a normalizar. Anteriormente, en el PCA no habíamos
normalizado, pues ya de por sí, escala los datos si activamos scale =
TRUE, y por tanto no afectaría a la visualización.

```{r normalizar}
# Guardar el conjunto train sin normalizar para la gráfica posterior
credit.Datos.Train.Original <- credit.Datos.Train

# Normalizar usando solo el conjunto de entrenamiento
preProc <- preProcess(credit.Datos.Train[, numeric_columns], method = c("center", "scale"))

# Aplicar la normalización al conjunto de entrenamiento
credit.Datos.Train[, numeric_columns] <- predict(preProc, 
                                        credit.Datos.Train[, numeric_columns])

# Aplicar la misma normalización al conjunto de prueba
credit.Datos.Test[, numeric_columns] <- predict(preProc, 
                                        credit.Datos.Test[, numeric_columns])

# Confirmar que los datos han sido normalizados
summary(credit.Datos.Train[, numeric_columns])
summary(credit.Datos.Test[, numeric_columns])

```

```{r}

# Gráfica de densidad antes de normalizar
d1 <- densityplot(
  formula(paste("~", paste(numeric_columns2, sep = "", collapse = " + "), collapse = "")),
  data = credit.Datos.Train.Original,
  main = "Variables sin Normalizar (train)",
  plot.points = FALSE
)

# Gráfica de densidad después de normalizar (usamos datos ya normalizados)
d2 <- densityplot(
  formula(paste("~", paste(numeric_columns2, sep = "", collapse = " + "), collapse = "")),
  data = credit.Datos.Train,
  main = "Variables Normalizadas (train)",
  plot.points = FALSE
)

# Mostrar ambas gráficas lado a lado
trellis.par.set(theme = standard.theme("pdf", color = FALSE))
print(d1, position = c(0, 0, 0.5, 1), more = TRUE)
print(d2, position = c(0.5, 0, 1, 1))
trellis.par.set(theme = standard.theme("pdf"))
```
Después de normalizar, todas las variables se transforman a un rango comparable. Esto se observa en el gráfico de variables normalizadas, donde las densidades están centradas alrededor de valores similares. Mientras que en el gráfico sin normalizar, vemos valores que llegan a decenas de miles, mientras que otros permanecen en rangos bajos

## Modelo rpart (Simple)


Rpart es un algoritmo basado en árboles de decisión en el que se va particionando de forma recursiva un conjunto de datos en subconjuntos de datos más homogéneos hasta que se cumpla un criterio de parada o no se obtengan mejores significativas. 

Para su correcto funcionamiento, se requiere de un conjunto de datos lo **suficientemente grande** y **representativo** del problema y se aconseja (aunque rpart pueda trabajar con nulos, outliers y sin escalar) un **preprocesamiento de los datos** previo (eliminar nulos, valores fuera de rango, normalizar/escalar...). Además:

* Realiza un control del tamaño y de la poda según el coste-complejidad del árbol, evitando así un sobreaprendizaje de los datos.

* Trabaja tanto con variables numéricas como con variables categóricas.

* Toma decisiones basándose en la **impureza**, cómo de concentrado está un subconjunto de datos en una sola clase. Escoge una división frente a otra si esta ofrece una menor impureza.

Uno de sus mayores puntos fuertes es su fácil interpretación y visualización de las relaciones entre las variables de entrada y la variable objetivo. Además, este modelo sirve como base para modelos más complejos y potentes como Random Forest.

Antes de empezar a entrenar el modelo, el método **train()** de caret con el que se entrenan los modelos **es capaz de encontrar la mejor configuración de hiperparametros** para nuestro problema y conjunto de datos. Además, con el método **trainControl()** de caret, podemos tener un mejor control sobre cómo va a realizar el entrenamiento (por ejemplo, para escoger la técnica de re-muestreo y configurarla). 

### Modelo con busqueda de hiperparámetros hecho por Caret

```{r}
# Preparar los datos
credit$V16 <- as.factor(credit$V16)
credit.Datos.Train$V16 <- as.factor(credit.Datos.Train$V16)
credit.Datos.Test$V16 <- as.factor(credit.Datos.Test$V16)


saveRpartCaret <- "ModeloRpartCaret.rds"

# Ajuste Manual mediante Grid Search
if (file.exists(saveRpartCaret)) {
  # Cargar el modelo manual si ya existe
  modelo_rpart_caret <- readRDS(saveRpartCaret)
  cat("Modelo Rpart cargado desde el archivo.\n")
}else{



# Configurar validación cruzada
set.seed(123)
control <- trainControl(
  method = "repeatedcv",  # Validación cruzada repetida
  number = 10,            # Número de particiones
  repeats = 3,            # Número de repeticiones
  verboseIter = FALSE      # Mostrar progreso
)


#Modelo con mejor configuración de hiperparámetros encontrada por caret.
set.seed(123)
modelo_rpart_caret <- train(
        V16 ~ ., 
        data = credit.Datos.Train,
        method = "rpart",
        trControl = control,
      )

# Guardar el modelo manual
  saveRDS(modelo_rpart_caret, saveRpartCaret)
  cat("Modelo Rpart entrenado y guardado.\n")

}

# Mostrar el arbol final
par(mar = c(1, 1, 1, 1))
rpart.plot(modelo_rpart_caret$finalModel,
           main = "Árbol de Decisión con Hiperparámetros encontrados por Caret", extra = 104)
par(mar = c(5, 4, 4, 2) + 0.1)


```

En cada nodo, nos muestran la clase predominante de los datos en ese nodo, la proporción de cada clase ( '-' y '+' respectivamente) y el porcentaje de datos en ese nodo respecto al conjunto de datos. **Vemos que la clase predominante en el conjunto de datos es '-'** (mostrado en el nodo raíz), además de que **V9** es la **variable más discriminante** en el conjunto de datos dado. 

Caret ha escogido la mejor combinación de hiperparámetros fijandose en la métrica **Acurracy**, que mide el número de predicciones correctas respecto al número de predicciones totales. Tambien muestra una métrica **Kappa**, que trata de representar lo mismo que Acurracy pero omitiendo los factores de azar, es decir, **trata de medir los aciertos que se provocan por cómo funciona el modelo**.

La combinación de hiperparámetros que Caret ha encontrado ha sido:

```{r}

##Imprimir combinaciones de hiperparametros que Caret ha probado, junto a su Acurracy y Kappa
print(modelo_rpart_caret) #Se queda con aquel con mejor Accuracy


```
### Modelo con busqueda de hiperparámetros hecho manualmente

Si nos fijamos en la documentación de rpart (<https://cran.r-project.org/web/packages/rpart/rpart.pdf>), Caret no trabaja con todos los hiperparámetros disponibles de rpart:

* **cp** (Complexity Parameter): Parámetro que indica la importancia relativa del tamaño del árbol con respecto al error. Cuanto mayor sea el valor de cp, mayor penalización por un mayor tamaño/compeljidad del árbol. **Parámetro usado para el control de la poda**. 

* **minsplit**: Número mínimo de ejemplos en un nodo para poder realizar su división. Valores más altos conllevan a que el árbol se quede con divisiones más generales de los datos.

* **minbucket**: Número mínimo de ejemplos en los nodos terminales. Altamente relacionada con el anterior hiperparámetro.

* **maxdepth**: Máxima profundidad del árbol. Valores más altos conllevan a una mayor posibilidad de sobreajuste si los anteriores hiperparámetros no han sido ajustados correctamente.

* Existen más parámetros, pero están relacionadas con la cross-validación o con el tratamiento en caso de valores nulos.

De todos estos, Caret **solo trabaja con Complexity Parameter**, sólo busca el mejor valor para este hiperparámetro. Por ello, **hemos decidido buscar manualmente la mejor configuración de hiperparametros**. Se exploran así combinaciones que están fuera del alcance de Caret y que pueden adaptarse mejor a nuestro problema y conjunto de datos.

* Para **cp**, se ha usado un rango de valores entre **0.001** y **0.01** por ser este el rango en el que caret ha buscado y obtenido mejores resultados.

* Para **minsplit**, se ha usado un rango de valores entre **10** y **30**, que, junto a Complexity Parameter nos permitirá buscar patrones de datos con un equilibrio entre especificas y generales. **minbucket se relaciona con minsplit** en cuanto a su valor, pues se recomiendan valores de minbuck iguales a la raíz cuadrada de minsplit y  un tercio de minsplit. Este rango de valores nos permite evitar divisiones especificas.

* Para **maxdepth**, se ha establecido un rango entre **5** y **15**. Valores inferiores a 5 podrían provocar un arbol muy general, mientras que 15 podrían ser muy específicos, presentando sobreajuste.

```{r}


#Modelo con mejor configuración de hiperparámetros encontrada de forma manual.


saveRpartManual <- "ModeloRpartManual.rds"

# Ajuste Manual mediante Grid Search
if (file.exists(saveRpartManual)) {
  # Cargar el modelo manual si ya existe
  modelo_y_resultados <- readRDS(saveRpartManual)
  results_df <- modelo_y_resultados$results_df
  results <- modelo_y_resultados$results
  modelo_rpart <- modelo_y_resultados$modelo
  cat("Modelo Rpart cargado desde el archivo.\n")
}else{

# Configurar validación cruzada
set.seed(123)
control <- trainControl(
  method = "repeatedcv",  # Validación cruzada repetida
  number = 10,            # Número de particiones
  repeats = 3,            # Número de repeticiones
  verboseIter = FALSE      # Mostrar progreso
)


# Inicializar para almacenar resultados
results <-list()

# Iterar sobre combinaciones de minsplit, minbucket, y maxdepth

for (minsplit in  c(10, 20, 30)) {
  for (minbucket in c(5, 7, 10)) {
    for (maxdepth in c(5, 10, 15)) {
      # Configurar control personalizado
      custom_control <- rpart.control(
        minsplit = minsplit,
        minbucket = minbucket,
        maxdepth = maxdepth
      )
      
      # Entrenar el modelo con train() y ajuste de cp
      set.seed(123)  # Para reproducibilidad
      modelo_rpart <- train(
        V16 ~ ., 
        data = credit.Datos.Train,
        method = "rpart",
        trControl = control,
        tuneGrid =  expand.grid(cp = seq(0.001, 0.1, by = 0.005)),
        control = custom_control  # Aplicar los valores personalizados
      )
      
      # Almacenar resultados
        results[[length(results) + 1]] <- list(
        cp = modelo_rpart$bestTune$cp,
        minsplit = minsplit,
        minbucket = minbucket,
        maxdepth = maxdepth,
        Accuracy = max(modelo_rpart$results$Accuracy),
        Kappa = max(modelo_rpart$results$Kappa),
        Modelo = modelo_rpart$finalModel
        )
      
    }
  }

}
# Convertir los resultados en un data.frame sin los modelos
results_df <- do.call(rbind, lapply(results, function(x) {
  data.frame(
    cp = x$cp,
    minsplit = x$minsplit,
    minbucket = x$minbucket,
    maxdepth = x$maxdepth,
    Accuracy = x$Accuracy,
    Kappa = x$Kappa
  )
}))


  # Guardar el modelo manual
  saveRDS(list(modelo = modelo_rpart, results_df=results_df, results=results), saveRpartManual)
  cat("Modelo RF entrenado y guardado.\n")


}


# Mostrar las combinaciones ordenadas por Accuracy
results_df <- results_df[order(-results_df$Accuracy), ]
head(results_df)

# Recuperar el mejor modelo
modelo_rpart_finalModel <- results[[which.max(results_df$Accuracy)]][["Modelo"]]

```

Podemos observar la diferencia entre Accuracy y Kappa entre los modelos cuya configuración de hiperparámetros ha sido buscada manualmente y con aquellas buscadas por Caret. **Como no conocemos nada sobre el conjunto de datos** , **no podemos determinar lo significativa que es esta diferencia** . Sin embargo, la mejor combinación encontrada manualmente presenta unos valores que provocan que este modelo sea más **simple** y **menos probable a presentar sobreajuste**.


```{r}
#Mostrar el arbol final
par(mar = c(1, 1, 1, 1))
rpart.plot(modelo_rpart_finalModel, main = 
          "Árbol de Decisión con Hiperparámetros manualmente",
          extra = 104)

par(mar = c(5, 4, 4, 2) + 0.1)

```


## Modelo rf (Complejo)

Random Forest es un algoritmo que se basa en árboles de decisión, en una combinación de **n** árboles de decisión **independientes**, obtenidos de la siguiente manera:

* Dado un conjunto de entrenamiento, se aplica **Bootstrapping** para generar **n** subconjuntos de entrenamiento (se recomienda como mínimo 100).
Cada uno de estos subconjuntos será utilizado por un único arból de decisión, dejando el resto de subconjuntos para validación.

* Sea **M** el número de variables de entrada del conjunto de entrenamiento original, se generr de forma aleatoria un número **m<< M**.

* Para cada árbol, se seleccionarán aletoriamente **m** variables diferentes, usadas para tomar decisiones a la hora de dividir un nodo. 

* Cada árbol predice una clase, siendo la **predicción final aquella clase** que **la mayoría de árboles han predicho**.

Al igual que rpart, para su correcto funcionamiento, se requiere de un conjunto de datos lo **suficientemente grande** y **representativo** del problema y se aconseja que se haya hecho previamente un **preprocesamiento de los datos**.

Con este funcionamiento se consigue entre otras cosas **evitar el sobreaprendizaje / overfitting **, obtener unas **predicciones con una alta precisión**. Sin embargo, su costo computacional es **alto** Atendiendo a la documentación de Caret (<https://topepo.github.io/caret/train-models-by-tag.html#random-forest>) sobre Random Forest, trabaja con el parámetro **mtry**, que se trata de **m**, es decir, el número de variables que cada árbol va a escoger de forma aleatoria.

### Modelo con busqueda de hiperparámetros hecha por Caret

Al tener en este conjunto de datos 15 variables, hemos considerado el rango de posibles valores de mtry para nuestro problema y conjunto de datos en [3,10].**Valores inferiores** pueden provocar que **no se detecten las relaciones entre las variables**, mientras que **valores superiores** pueden provocar **sobreajuste**, por lo que hemos considerado innecesaria la busqueda por todo el rango posible de valores.

```{r}
#Busuqeda de hiperparámetros por Caret

saveRFCaret <- "ModeloRFCaret.rds"

# Ajuste Manual mediante Grid Search
if (file.exists(saveRFCaret)) {
  # Cargar el modelo manual si ya existe
  modelo_rf_caret <- readRDS(saveRFCaret)
  cat("Modelo RF cargado desde el archivo.\n")
}else{


# Configurar validación cruzada
set.seed(123)
control <- trainControl(
  method = "repeatedcv",  # Validación cruzada repetida
  number = 10,            # Número de particiones
  repeats = 3,            # Número de repeticiones
  verboseIter = FALSE      # Mostrar progreso
)

# Entrenar el modelo Random Forest
set.seed(123)
modelo_rf_caret <- train(
  V16 ~ ., 
  data = credit.Datos.Train, 
  tuneGrid = expand.grid(mtry = c(3,5,7,10)),
  method = , 
  trControl = control, 
)

# Guardar el modelo manual
  saveRDS(modelo_rf_caret, saveRFCaret)
  cat("Modelo RF entrenado y guardado.\n")

}
# Resultados del modelo
print(modelo_rf_caret)

#La mejor combinacion de hiperparametros encontrado por caret ha sido el siguiente:
print(modelo_rf_caret$bestTune)


```

### Modelo con busqueda de hiperparámetros hecho manualmente

Según la documentación de Random Forest en R (<https://cran.r-project.org/web/packages/randomForest/randomForest.pdf>) existen varios parámetros más, entre el que se destaca **ntree** y **nodesize**, es decir, el número de árboles que tendrá el bosque (por defecto, es igual a 500) y el tamaño mínimo de los nodos terminales (por defecto, es igual a 1).

* **ntree** es muy influyente, puesto que, a mayor coste, **mayor** va a ser su **coste computacional**, pero **mejor precisión** va a obtener (ojo, tamaños excesivos pueden provocar aún asi sobreajuste).

* **nodesize** tambien influye bastante, puesto que a mayor número mínimo, mayor la probabilidad de que el árbol aprenda patrones generales de los datos, evitando el sobreajuste.


```{r}

saveRFManual <- "ModeloRFManual.rds"

if (file.exists(saveRFManual)) {
  # Cargar el modelo manual si ya existe
  modelo_y_resultados <- readRDS(saveRFManual)
  results_df <- modelo_y_resultados$results_df
  results <- modelo_y_resultados$results
  modelo_rf <- modelo_y_resultados$modelo
  
  cat("Modelo manual RF cargado desde el archivo.\n")
}else{


# Configurar validación cruzada
set.seed(123)
control <- trainControl(
  method = "repeatedcv",  # Validación cruzada repetida
  number = 10,            # Número de particiones
  repeats = 3,            # Número de repeticiones
  verboseIter = FALSE      # Mostrar progreso
)

#Busqueda manual de hiperparámetros
mtrys <- c(3, 5, 7, 10)
nodesizes <- c(5, 7, 10)
ntrees <- c(500,1000, 1500)
  
  results <-list()
  for(nodesize in nodesizes){
   for (ntree in ntrees) {
  
      set.seed(123)
      modelo_rf <- train(
          V16 ~ .,
          data = credit.Datos.Train,
          method="rf",
          tuneGrid = expand.grid(mtry = mtrys),
          trControl=control,
          ntree=ntree,
          nodesize = nodesize)
    
       results[[length(results) + 1]] <- list(
          mtry = modelo_rf$bestTune$mtry,
          ntree = ntree,
          nodesize = nodesize,
          Accuracy = max(modelo_rf$results$Accuracy),
          Kappa = max(modelo_rf$results$Kappa),
          Modelo = modelo_rf$finalModel
        )

  
   }
  }
       
        # Convertir los resultados en un data.frame sin los modelos
      results_df <- do.call(rbind, lapply(results, function(x) {
      data.frame(
       mtry = x$mtry,
       ntree = x$ntree,
       nodesize = x$nodesize,
      Accuracy = x$Accuracy,
      Kappa = x$Kappa
      )
      }
      )
      )
  
# Guardar el modelo manual
  saveRDS(list(modelo = modelo_rf, results_df=results_df, results=results), saveRFManual)
  cat("Modelo RF entrenado y guardado.\n")

}


  results_df <- results_df[order(-results_df$Accuracy), ]
  modelo_rf_finalModel <-  results[which.max(
  results_df$Accuracy)]["Modelo"]
  
  
  print(results_df)


```


Al comparar el modelo generado automáticamente por Caret y los obtenidos mediante la búsqueda manual de hiperparámetros, observamos que las combinaciones de hiperparámetros encontradas de forma manual obtienen un mejor rendimiento respecto a los que encuentra Caret. De nuevo, **como no tenemos contexto sobre el conjunto de datos**, **no podemos saber si esta diferencia se considera significativa o no**. Es por eso que se elige la configuración **mtry = 7 ** y **ntree = 500** y **nodesize = 7** , al tener valores que provocan que los árboles generados sean **más generales** y tengan un **menor coste computacional**.




## Modelo nnet (Complejo)

Este modelo trata de imitar el funcionamiento de las redes neuronales del cerebro humano. **Se compone de una serie de capas**, **conectadas entre sí**:

-   Una **capa de entrada** que contiene aquellas **variables de entrada
    predictoras**.

-   Una **capa oculta** encargada de procesar los datos pasados por la **capa
    de entrada**. Esto lo realiza mediante

-   Una **capa de salida** en el que se genera el **resultado final de una
    predicción**.

Como dice la documentación
(<https://rubenfcasal.github.io/aprendizaje_estadistico/implementaci%C3%B3n-en-r-2.html>),
el modelo nnet implementa redes neuronales de tipo feedforward con una única capa oculta. Este tipo de red neuronal es **útil para resolver problemas de clasificación y regresión** gracias a su **capacidad para modelar relaciones no lineales complejas** entre las variables de entrada y la salida. Si queremos implementar redes neuronales más complejas con múltiples capas, debemos usar paquetes como keras o tensorflow.

Antes de empezar, es fundamental **escalar** las **variables predictoras** para que estén en **rangos similares**, ya que la red neuronal es muy sensible y los pesos iniciales se optimizan mejor cuando las variables tienen rangos comparables. También es necesario que las variables **categóricas** se conviertan en **variables dummy**. Anteriormente, en el PCA hicimosvariables dummy, pero las vamos a volver hacer para tenerlo por separado, y que sean específicamente para entrenar el modelo.

```{r}
objetivo_train <- credit.Datos.Train$V16
objetivo_test <- credit.Datos.Test$V16

predictoras_train <- credit.Datos.Train[, names(credit.Datos.Train) != "V16"]
predictoras_test <- credit.Datos.Test[, names(credit.Datos.Test) != "V16"]

categorical_columns <- names(predictoras_train)[sapply(predictoras_train, is.factor)]

# Crear dummies con modelo de entrenamiento
dummies_model <- dummyVars(~ ., data = predictoras_train[, categorical_columns, drop = FALSE], fullRank = TRUE)
dummies_train <- predict(dummies_model, newdata = predictoras_train[, categorical_columns, drop = FALSE])
dummies_test <- predict(dummies_model, newdata = predictoras_test[, categorical_columns, drop = FALSE])

predictoras_train_dummies <- cbind(
  predictoras_train[, sapply(predictoras_train, is.numeric), drop = FALSE],
  dummies_train
)
predictoras_test_dummies <- cbind(
  predictoras_test[, sapply(predictoras_test, is.numeric), drop = FALSE],
  dummies_test
)

# Estas ya están normalizadas (porque aplicamos preProc antes)
# Añadir V16
credit_dummies_train_normalizado <- cbind(predictoras_train_dummies, V16 = objetivo_train)
credit_dummies_test_normalizado <- cbind(predictoras_test_dummies, V16 = objetivo_test)

# Asegurarnos que V16 es factor con mismos niveles
credit_dummies_train_normalizado$V16 <- factor(credit_dummies_train_normalizado$V16, levels=c("+","-"))
credit_dummies_test_normalizado$V16 <- factor(credit_dummies_test_normalizado$V16, levels=c("+","-"))


```

Por otro lado, el modelo también se ve **favorecido** con la **reducción de dimensionalidad**, ya que pueden sobreajustarse o perder eficiencia si los datos tienen muchas variables irrelevantes o correlacionadas. Por lo que usaremos el PCA calculado anteriormente.

Viendo la documentación de Caret, para redes neuronales
(<https://topepo.github.io/caret/train-models-by-tag.html#neural-network>),
recopilamos los hiper-parámetros más importantes que afectan el
rendimiento del modelo:

-   **size (Número de nodos en la capa oculta):** Controla la capacidad
    del modelo para aprender patrones complejos. Valores bajos pueden
    subestimar los patrones de los datos, mientras que valores altos
    pueden provocar sobreajuste.

-   **decay (Regularización L2):** Penaliza los pesos altos en la red,
    lo que ayuda a prevenir el sobreajuste. Valores pequeños (como
    0.001) suelen ser adecuados; sin embargo, pueden ser necesarios
    valores mayores para datos ruidosos.

-   **maxit (Iteraciones máximas):** Controla el número máximo de
    iteraciones del algoritmo de optimización. Por defecto, es 100, pero
    puede aumentarse si el modelo no converge.

-   **linout (Salida lineal):** Se establece en TRUE para problemas de
    regresión y en FALSE (por defecto) para clasificación.

-   **MaxNWts (Número máximo de pesos):** Define el número total de
    pesos que la red puede manejar. Por defecto es 1000, pero debe
    aumentarse si el número de variables o nodos es alto.

-   **Funciones de activación:** Por defecto, nnet utiliza una función
    sigmoide para problemas de clasificación binaria.


### Ajuste automático de hiper-parámetros con Caret

Las redes neuronales son modelos **altamente sensibles a los hiperparámetros**, y una **mala selección inicial** de estos puede comprometer **significativamente** el rendimiento del modelo. Por ello, **hemos decidido comenzar con un ajuste automático**. Caret prueba automáticamente 10 combinaciones de los hiperparámetros(size y decay), usando la cantidad especificada en tuneLength. Para cada
combinación, realiza validación cruzada con las particiones definidas en trainControl, en este caso 10. Después, elige automáticamente la combinación con el mejor rendimiento según Accuracy y Kappa.

```{r}

# Definir control para validación cruzada
control <- trainControl(
  method = "cv",            # Método de validación cruzada
  number = 10,              # Número de particiones (folds)
  verboseIter = FALSE        # Mostrar progreso
)

# Verificar y entrenar el modelo automático con Caret
if (!file.exists("model_nnetCaret.rds")) {
# Entrenar el modelo de red neuronal con ajuste automático
  set.seed(123)
  model_nnetCaret <- train(
    V16 ~ .,                   # Fórmula: variable objetivo frente a todas las variables predictoras
    data = credit_dummies_train_normalizado, # Datos normalizados con dummies
    method = "nnet",           # Modelo de red neuronal
    trControl = control,       # Validación cruzada
    tuneLength = 10,           # Ajustar automáticamente 10 combinaciones de hiperparámetros
    linout = FALSE,            # Usar función logística para clasificación binaria
    trace = FALSE,             # No mostrar mensajes de convergencia
    MaxNWts = 10000            # Límite de pesos en la red (para muchas variables)
  )
  saveRDS(model_nnetCaret, "model_nnetCaret.rds")
} else {
  model_nnetCaret <- readRDS("model_nnetCaret.rds")
}

# Imprimir resultados del modelo
print(model_nnetCaret)

# Ver la mejor combinación de hiperparámetros
cat("Mejores hiperparámetros encontrados:\n")
print(model_nnetCaret$bestTune)

```
En este caso, los mejores hiperparámetros son **size=1** y **decay=0.007498942**
Dando un **accuracy** de **0.8659412** y un **Kappa** de **0.7248953**

### Búsqueda manual de hiper-parámetros

En este apartado, realizaremos un ajuste manual de los hiperparámetros 
del modelo, basándonos en los valores obtenidos durante el ajuste automático 
previo. Nuestro objetivo es refinar los resultados iniciales, explorando 
un rango más acotado de combinaciones que son favorables, con tal de mejorar 
el rendimiento del modelo ,y encontrar configuraciones que maximicen el 
Accuracy y el Kappa.

Nos centraremos en **dos hiperparámetros clave**:

- **Número de nodos en la capa oculta** (**size**): Vemos que el valor 1 es el más 
prometedor y con el mayor Accuracy, y que con valores mayor a 3 tiende a 
reducir un poco el rendimiento. Por lo que probaremos con valores entre 1 y 3.

- **Tasa de regularización** (**decay**): El mejor valor está cerca de valores como 
0.0075 Por lo que probaremos con valores cercanos a estos.


```{r}

# Configuración de la validación cruzada
control <- trainControl(
  method = "cv",           # Validación cruzada
  number = 10,             # 10 pliegues
  verboseIter = FALSE       # Mostrar progreso
)


if (!file.exists("nnet_modelManual.rds")) {
  # Configuración del grid de hiperparámetros
  nnet_grid <- expand.grid(
      size = c(1, 2, 3),                     # Focalizado en los valores más prometedores
      decay = decays <- c(0.00000, 0.00001,0.00010, 0.000237,  # Valores intermedios entre potencias de 10
              0.000562,0.00100,0.00200,0.00500,0.01000,0.02000,0.05000,0.10000)
  )

  # Entrenar el modelo usando train() y tuneGrid
  set.seed(123)
  nnet_modelManual <- train(
    V16 ~ .,                       # Fórmula del modelo
    data = credit_dummies_train_normalizado, # Dataset con dummies normalizado
    method = "nnet",               # Modelo: Red neuronal
    trControl = control,           # Validación cruzada
    tuneGrid = nnet_grid,          # Grid de hiperparámetros
    linout = FALSE,                # Clasificación binaria (no salida lineal)
    trace = FALSE,                 # No mostrar mensajes internos de optimización
    MaxNWts = 10000                # Límite de pesos (para redes grandes)
  )
  saveRDS(nnet_modelManual, "nnet_modelManual.rds")
} else {
  nnet_modelManual <- readRDS("nnet_modelManual.rds")
}

# Mostrar el modelo entrenado
print(nnet_modelManual)

# Ver la mejor combinación de hiperparámetros
cat("Mejores hiperparámetros encontrados:\n")
print(nnet_modelManual$bestTune)

```

De esta manera, obtenemos que **los mejores resultados** son con un **size=1** y un **decay=1e-04**	Dando un **Accuracy** de **0.8795158** y un **Kappa** de **0.7530687**.

### Ajuste automático hiper-parámetros con PCA

En este apartado, implementaremos el PCA calculado anteriormente. Nosotros vamos 
volver a normalizar las variables, a pesar de que el PCA las genera escaladas. Pero,
queremos garantizar que estén completamente normalizadas, y asegurarnos de que
no haya problemas para el modelo, pues como mencionamos anteriormente, es importante
este requisito. 

Por otra parte, **al aplicar otro preprocesado**, **no nos sirven los hiperparámetros del grid anterior**. Las características han cambiado y **al aplicar PCA**, **las variables ya no son las mismas**, su escala y correlaciones cambian, y esto puede alterar 
las necesidades del modelo. Por ello, debemos realizar un **ajuste independiente** de hiperparámetros **tras aplicar PCA**.


Vamos a volver a probar, a realizar un ajuste automático de los hiper-parámetros size y decay. Posteriormente, volveremos a realizar un ajuste manual con los resultados obtenidos, con tal de ver si conseguimos sacar un mayor acuraccy.

```{r}

# Normalizar
credit.Datos.Train.Transf.PCA <- scale(credit.Datos.Train.Transf.PCA)
credit.Datos.Test.Transf.PCA <- scale(credit.Datos.Test.Transf.PCA)

# Datos tras el PCA
train_pca <- credit.Datos.Train.Transf.PCA
test_pca <- credit.Datos.Test.Transf.PCA

# Objetivo como factor
train_target <- as.factor(credit.Datos.Train$V16)
test_target <- as.factor(credit.Datos.Test$V16)

# Configurar validación cruzada
control <- trainControl(
  method = "cv",         # Validación cruzada
  number = 10,           # Número de particiones
  verboseIter = FALSE     # Mostrar progreso
)

# Verificar y entrenar modelo automático con PCA
if (!file.exists("model_nnet_auto_pca.rds")) {
  set.seed(123)
  model_nnet_auto_pca <- train(
    x = train_pca,
    y = train_target,
    method = "nnet",
    trControl = control,
    tuneLength = 10,        # Ajuste automático con 10 combinaciones
    linout = FALSE,
    trace = FALSE,
    MaxNWts = 10000
  )
  saveRDS(model_nnet_auto_pca, "model_nnet_auto_pca.rds")
} else {
  model_nnet_auto_pca <- readRDS("model_nnet_auto_pca.rds")
}




print(model_nnet_auto_pca)
cat("Mejores hiperparámetros (con PCA, automático):\n")
print(model_nnet_auto_pca$bestTune)

```

Bajo el escenario con PCA, la **mejor combinación** hallada durante el ajuste automático fue **size** = **1** y **decay** = **0.04216965**, con un **Acuraccy** de **0.8387119** y un **Kappa** de **0.6689545**

### Búsqueda manual de hiper-parámetros con PCA

Seguidamente utilizaremos un grid refinado con los resultados del modelo anterior, que juntaremos también con el grid del modelo sin PCA, para realizar una búsqueda manual. De esta manera, **con el PCA vamos a evaluar si la reducción de dimensionalidad mejora el rendimiento del modelo** . Esto también nos permitirá explorar con mayor precisión.

```{r}

if (!file.exists("model_nnet_manual_pca.rds")) {
  nnet_grid_pca <- expand.grid(
    size = c(1, 2, 3, 5),
    decay = decays <- c(0.000000, 0.000001, 0.000005, 0.000010, 0.000020, 0.000050, 0.000100, 0.000200, 0.000237, 
              0.000300, 0.000500, 0.000562, 0.001000, 0.002000, 0.004000, 0.004500, 0.005000, 0.005500, 
              0.006000, 0.006500, 0.006800, 0.007000, 0.007100, 0.007200, 0.007300, 0.007400, 0.007500, 
              0.007600, 0.007700, 0.007800, 0.007900, 0.008000, 0.008500, 0.009000, 0.009500, 0.010000, 
              0.020000, 0.050000, 0.100000)
  )
  
  set.seed(123)
  model_nnet_manual_pca <- train(
    x = train_pca,
    y = train_target,
    method = "nnet",
    trControl = control,
    tuneGrid = nnet_grid_pca,
    linout = FALSE,
    trace = FALSE,
    MaxNWts = 10000
  )
  saveRDS(model_nnet_manual_pca, "model_nnet_manual_pca.rds")
} else {
  model_nnet_manual_pca <- readRDS("model_nnet_manual_pca.rds")
}

print(model_nnet_manual_pca)
cat("Mejores hiperparámetros (con PCA, manual):\n")
print(model_nnet_manual_pca$bestTune)


```

Sorprendentemente, hemos conseguido mejorar un poco más con un **size**=**1** y un 
**decay**=**0.01**, dando un nuevo **Accuracy** de **0.8467843** y un **Kappa de **0.6877926**.

### Comparativa final

El primer modelo, sin reducción de dimensionalidad vemos que tiene un accuracy
mayor al obtenido con el PCA. No obstante, vamosa compararlos uno a uno en cuanto a tiempo de ejecución, con tal de considerar si merece la pena la reducción o no.

```{r}
# Mejores hiperparámetros sin PCA
mejor_grid_no_pca <- nnet_modelManual$bestTune 

# Mejores hiperparámetros con PCA
mejor_grid_pca <- model_nnet_manual_pca$bestTune

# Tiempo de ejecución sin PCA
start_time_no_pca <- Sys.time()
set.seed(123)
nnet_model_no_pca <- train(
  V16 ~ .,
  data = credit_dummies_train_normalizado,
  method = "nnet",
  trControl = control,
  tuneGrid = mejor_grid_no_pca,
  linout = FALSE,
  trace = FALSE,
  MaxNWts = 10000
)
end_time_no_pca <- Sys.time()
execution_time_no_pca <- end_time_no_pca - start_time_no_pca

# Tiempo de ejecución con PCA
start_time_pca <- Sys.time()
set.seed(123)
nnet_model_pca <- train(
  x = train_pca,
  y = train_target,
  method = "nnet",
  trControl = control,
  tuneGrid = mejor_grid_pca,
  linout = FALSE,
  trace = FALSE,
  MaxNWts = 10000
)
end_time_pca <- Sys.time()
execution_time_pca <- end_time_pca - start_time_pca

cat("Tiempo de ejecución SIN PCA:", execution_time_no_pca, "\n")
cat("Tiempo de ejecución CON PCA:", execution_time_pca, "\n")
```

Aunque el **modelo con PCA** es **ligeramente más rápido**,la pérdida en Accuracy de **0.8795158** a **0.8467843** es **bastante notoria**.

Es verdad que la reducción de dimensionalidad beneficia al desempeño del modelo, pero nosotros **consideramos que no merece la pena esa pérdida en el Accuracy** por tan solo unas milésimas de segundo. Como sabemos, el **database** que estamos empleando ya está **bastante trabajado**, en el sentido de que **si quitamos variables**, está claro que **va a bajar el accuracy** y que realmente todas las variables tienen su importancia a la hora de predecir. Si este fuese más extenso, con más ruido y con más variables irrelevantes (realmente hay poca correlación),  el PCA podría haber ayudado más. Aquí, no merece la pena la reducción de dimensionalidad, pues no se ha observado una ventaja clara en el rendimiento ni en el tiempo. Por lo tanto, en cuanto a redes neuronales, vamos a escoger el modelo sin PCA **nnet_modelManual** , que emplea como hiperparámetros **size**=**1**, **decay**=**1e-04**, y logrando un **Acuraccy** de **0.8795158**.


## Modelo gbm (Complejo)

**GBM funciona combinando modelos débiles**, **como los árboles de decisión**, genera secuencialmente árboles que van corrigiendo los errorees de los anterioeres. Así **mejora de forma progresiva** el resultado **hasta llegar a su predicción final**.

```{r}
modelLookup("gbm")

```

Como podemos observar, todos los parámetos son aptos tanto para
regresión como clasificación y también puede generar probabilidades
asociadas a las clases

Consultando con la documentación en R de gbm
(<https://search.r-project.org/CRAN/refmans/traineR/html/train.gbm.html>)
encontramos los siguientes hiperparámetros:

-   **n.tree**s: define el número de modelos debiles que se crearán. Cuanto
    más esto es el valor, más capaz es de aprender patrones, pero
    también aumenta mucho el coste computacional además de se más
    propenso al overfitting.

-   **interaction.deph**: es la profundidad máxima que puede alcanzar un
    árbol generado. Define la complejidad de los árboles generados y su
    habilidad para indentificar interacciones entre variables
    predictoras

-   **shrinkage**: es la tasa de aprendizaje del modelo. Pondera la
    contribución de cada árbol sobre el resultado final. Los valores
    bajos suelen dar un resultado más preciso, aunque precisan de más
    árboles.

-   **n.minobsinnode**: es el número mínimo de datos en un nodo terminal del
    árbol. En este caso, los valores más altos suelen dar una mejor
    generalización.

A diferencia de los dos primeros modelos, en el caso de que se tengan
**variables categóricas**, **se requieren transformar a variables dummy*.**
Además, **gbm** es un modelo que tiene una **alta demanda computacional**, sobre
todo con conjuntos de datos grandes o cuando ajustamos muchos parámetros. También es un modelo que tiende a sufrir de **overfitting si no se ajustan bien los parámetros**.

```{r}
credit.Datos.Train$V16 <- as.factor(credit.Datos.Train$V16)
credit.Datos.Test$V16 <- as.factor(credit.Datos.Test$V16)
# Configuración de los archivos donde se guardarán los modelos
model_auto_path <- "model_gbm_auto.rds"
model_grid_path <- "model_gbm_grid.rds"
model_grid_fine_path <- "model_gbm_fine_tuned.rds"
```

### Ajuste Automático

**Caret prueba con un grid por defecto automáticamente 10 combinaciones** (se lo indicamos con tuneLength) de los hiperparámetros de **interaction.depth** y **n.trees**, manteniendo **shrinkage** y **n.minobsinnode constante**, usando la cantidad especificada en tuneLength. Para cada combinación, realiza validación cruzada con las particiones definidas en trainControl, en este caso 10, y a 3 repeticiones. Después, elige automáticamente la combinación con el mejor rendimiento **según Accuracy y Kappa**.

```{r}
# Ajuste Automático
if (file.exists(model_auto_path)) {
  # Cargamos el modelo automático si ya existe
  model_gbm_auto <- readRDS(model_auto_path)
  cat("Modelo automático cargado desde el archivo.\n")
} else {
  set.seed(123)
  control_auto <- trainControl(
    method = "repeatedcv",    # Cross validation
    number = 10,              # Dividir los datos en 10 folds
    repeats = 3,              # Repetir 3 veces
    verboseIter = FALSE,      # No mostrar progreso durante el ajuste
  )
  
  # Entrenamos el modelo GBM con el conjunto de datos transformado y normalizado
  set.seed(123)
  model_gbm_auto <- train(
    V16 ~ .,                   # Fórmula: variable objetivo frente a todas las predictoras
    data = credit.Datos.Train, # Usar el conjunto de datos de entrenamiento
    method = "gbm",            # Método GBM
    trControl = control_auto,  # Control de entrenamiento personalizado
    tuneLength = 10,           # Probar 10 combinaciones aleatorias de hiperparámetros
    verbose = FALSE            # Desactivar mensajes de progreso
  )

 # Guardamos el modelo
  saveRDS(model_gbm_auto, model_auto_path)
  cat("Modelo automático entrenado y guardado.\n")
}

# Imprimimos los resultados del modelo automático (ver .Rmd)
# print(model_gbm_auto)

best_gbm_auto <- model_gbm_auto$bestTune
cat("\nMejores hiperparámetors modelo GBM automático:\n")
print(best_gbm_auto)

cat("\nAccuracy mejor modelo GBM automático entrenado:\n")
print(max(model_gbm_auto$result$Accuracy))

```

Vemos que la mejor configuración conseguida es **n.trees**=**50**, **interaction.depth**=**5**, **shrinkage**=**0.1**, **n.minobsinnode**=**10**.

### Ajuste manual con grid

Ahora haremos un ajuste manual con grid. Ajustaremos las mismas
variables pero esta vez definiendo **nosotros** el rango de búsqueda de los
hiperparámetros con grid. Indicamos en el control que quedemos hacer una
busqueda con grid, creamos el grid con **expand.grid()** y los parámetros que
nosotros decidimos, y se lo damos a **train()** en el parámetro **tuneGrid**.

```{r}
# Ajuste Manual mediante Grid Search
if (file.exists(model_grid_path)) {
  # Cargamos el modelo manual si ya existe
  model_gbm_grid <- readRDS(model_grid_path)
  cat("Modelo manual con Grid cargado desde el archivo.\n")
} else {
  # Entrenamos el modelo manual si no existe
  set.seed(123)
  control_manual <- trainControl(
    method = "repeatedcv",       # Cross validation
    number = 10,                 # Dividir los datos en 10 folds
    repeats = 3,                 # Repetir 3 veces
    verboseIter = FALSE,         # No mostrar progreso
    search = "grid"              # Ajuste de hiperparámetros con grid
  )
  
  # Definimos un grid de hiperparámetros personalizado
  grid_manual <- expand.grid(
    n.trees = c(30, 50, 100, 200, 500),             # Números de árboles a probar
    interaction.depth = c(3, 5, 8, 10, 12),         # Profundidades de árbol a probar
    shrinkage = c(0.01, 0.05, 0.1, 0.15, 0.2),      # Tasas de aprendizaje a probar
    n.minobsinnode = c(5, 10, 15, 20)               # Número mínimo de observaciones en un nodo a probar
  )
  
  
# Entrenamos el modelo con grid que acababmos de definir
  set.seed(123)
  model_gbm_grid <- train(
    V16 ~ .,                    # Fórmula: variable objetivo frente a todas las predictoras
    data = credit.Datos.Train,  # Usar el conjunto de datos de entrenamiento
    method = "gbm",             # Método GBM
    trControl = control_manual, # Control de entrenamiento
    tuneGrid = grid_manual,     # Hiperparámetros definidos en el grid
    verbose = FALSE             # Desactivar mensajes de progreso
  )
# Guardamos el modelo manual
  saveRDS(model_gbm_grid, model_grid_path)
  cat("Modelo manual (Grid Search) entrenado y guardado.\n")
}
# Imprimimos los resultados del modelo manual (ver .Rmd)
# print(model_gbm_grid)

best_gbm_grid <- model_gbm_grid$bestTune
cat("\nMejores hiperparámetors modelo GBM manual:\n")
print(best_gbm_grid)

cat("\nAccuracy mejor modelo GBM manual entrenado:\n")
print(max(model_gbm_grid$result$Accuracy))
```

Con el ajuste de grid hemos hecho una búsqueda grusa, y posteriormente afinaremos la búsqueda en baso a los resultados. En este caso hemos obtenido la configuración **n.trees**=**50**, **interaction.depth**=**8**, **shrinkage**=**0.1**, **n.minobsinnode**=**10**.

### Ver la mejor configuración que ha encontrado cada uno

**1.- Boxplots de Accuracy**

```{r}
# Preparamos los datos para el gráfico
results_auto <- model_gbm_auto$results
results_grid <- model_gbm_grid$results
results_auto$Method <- "Automático"
results_grid$Method <- "Grid personalizado"

combined_results <- bind_rows(results_auto, results_grid)

# Crear boxplot
ggplot(combined_results, aes(x = Method, y = Accuracy, fill = Method)) +
  geom_boxplot() +
  stat_summary(fun = mean, geom = "text", aes(label = round(after_stat(y), 3)), vjust = -0.5) +
  theme_minimal() +
  labs(
    title = "Comparación de accuracy entre Métodos",
    x = "Método",
    y = "Accuracy",
    fill = "Método"
  )
```

La configuración con grid ha obtenido un promedio de **Accuracy** **un poco mejor** (**0.873**) **comparándolo con el método de ajuste automático** (**0.871**). Sin embargo, la diferencia es realmente pequeña, podemos decir que que ambos métodos tienen de media tienen un rendimiento similar.

En cuanto a la variabilidad, son muy parecidos en ambos casos, aunque el automático es ligeramente menor. Esto se aprecia por su rango intercuartil más estrecho. Por contra, el método grid presenta ligeramente más dispersión en los resultados, lo que quiere decir que explora más configuraciones y diferentes.

Además vemos que el método grid, a diferencia del automático, sí genera algunos valores atípicos. Estos pueden ser combinaciones específicas que no se han ajustado bien.


**2.- Curvas de Desempeño Promedio**

Una gráfica que hemos considerado interesante es la del **desempeño promedio**, que nos **permite analizar cómo afecta cada hiperparámetro** al **accuracy** del modelo e identifica los valores óptimos de los mismos y algunas de las tendencias al ajustarlos como el **overfitting**. Además, nos **ayuda a comparar los diferentes ajustes de hiperparámetros** probados en términos de **estabilidad** y capacidad de exploración del espacio de **hiperparámetros**. También nos es muy útil para guiar el **Fine Tuning** en rangos específicos que realizaremos después.


```{r}
# Listamos los hiperparámetros testeados
variables <- c("n.trees", "interaction.depth", "shrinkage", "n.minobsinnode")

# Creamos los gráficos para cada parámetro
plots <- lapply(variables, function(var) {
  
  # Agrupamos los datos por hiperparámetors
  auto_grouped <- group_by(results_auto, across(all_of(var)))
  # Caluculamos las medias de cada grupo y vamos limpiando los grupos
  auto_summarized <- summarize(auto_grouped, Accuracy = mean(Accuracy), .groups = 'drop')
  # Agregamos la columna modelo para poder indentificarlos luego
  auto_avg <- mutate(auto_summarized, Modelo = "Automático")

  # Repetimos lo anterior para el modelo con grid manual
  grid_grouped <- group_by(results_grid, across(all_of(var)))
  grid_summarized <- summarize(grid_grouped, Accuracy = mean(Accuracy), .groups = 'drop')
  grid_avg <- mutate(grid_summarized, Modelo = "Grid personalizado")
  
  # Combinamos los resultados obtenidos
  combined_avg <- bind_rows(auto_avg, grid_avg)
  
  # Creamos el gráfico final
  ggplot(combined_avg, aes_string(x = var, y = "Accuracy", color = "Modelo", group = "Modelo")) +
    geom_line(size = 1) +
    geom_point(size = 2) +
    labs(
      title = paste("Impacto de", var, "sobre Accuracy"),
      x = var,
      y = "Accuracy",
      color = "Modelo"
    ) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)
  )
})

# Mostramos los gráficos en conjunto
grid.arrange(grobs = plots, ncol = 2)


```

Podemos ver que **cuanto mayor** es el valor de **n.trees**, **más bajo** es el promedio de **Accuracy**. Esto quiere decir que posiblemente un valor alto en este parámetro lleva al **sobreajuste**. Esto lo observamos también para el resto de parámetros. También observamos que el método **grid** es más consistente gracias a que **explora muchas posibilidades**, mientras que el automático es más caótico. Además podemos destacar que el método automático ajusta solo **n.trees** e **interaction.depth**. Ambas configuraciones muestran su pico más alto alrededor del valor **50** en **n.trees** y picos buenos entre **5** y **10** para **interaction.depth**, lo que nos indica que son buenos valores para entrenar modelos puesto que están en el pico del rendimiento antes de sobreajustar. En cuanto a **shrinkage** vemos que el valor óptimo sería alrededor de **0.05** y para **n.minobsinnode** sería cerca de **10**.

### Fine Tune

Vamos a hacer a continuación un fine tunning con las conclusiones que hemos obtenido para ver si conseguimos mejorar el rendimiento de nuestro modelo.

```{r}
# Verificamos si el modelo ya existe en disco
if (file.exists(model_grid_fine_path)) {
  # Cargamos el modelo desde el archivo RDS si existe
  model_gbm_fine <- readRDS(model_grid_fine_path)
  cat("Modelo cargado desde archivo:", model_grid_fine_path, "\n")
} else {
  # Definimos el grid de hiperparámetros con un ajuste más fino
  grid_fine_tuning <- expand.grid(
    n.trees = seq(40, 60, by = 5),             # Alrededor de 50
    interaction.depth = seq(5, 10, by = 1),    # Entre 5 y 10
    shrinkage = seq(0.03, 0.1, by = 0.01),    # Alrededor de 0.05
    n.minobsinnode = seq(8, 12, by = 1)        # Valores cercanos a 10
  )

  set.seed(123)
  control_fine <- trainControl(
    method = "repeatedcv",       # Cross validation
    number = 10,                 # Dividir los datos en 10 folds
    repeats = 3,                 # Repetir 3 veces
    verboseIter = FALSE,         # No mostrar progreso
    search = "grid"              # Ajuste de hiperparámetros con grid
  )

  # Entrenamos el modelo GBM con el ajuste fino
  set.seed(123)
  model_gbm_fine <- train(
    V16 ~ .,                      # Fórmula: variable objetivo frente a todas las predictoras
    data = credit.Datos.Train,    # Usar el conjunto de datos de entrenamiento
    method = "gbm",               # Método GBM
    trControl = control_fine,     # Control de entrenamiento
    tuneGrid = grid_fine_tuning,  # Grid de hiperparámetros con ajuste fino
    verbose = FALSE               # Desactivar mensajes de progreso
  )

  # Guardar el modelo en formato RDS
  saveRDS(model_gbm_fine, model_grid_fine_path)
  cat("Modelo guardado en archivo:", model_grid_fine_path, "\n")
}

# Imprimimos los resultados del modelo manual (ver .Rmd)
# print(model_gbm_fine)

# Imprimir el mejor modelo encontrado
cat("\nMejores hiperparámetors modelo GBM fine tune:\n")
print(model_gbm_fine$bestTune)

cat("\nAccuracy mejor modelo GBM fine tune entrenado:\n")
print(max(model_gbm_fine$result$Accuracy))
```

La mejor configuración encontrada en este caso es **n.trees**=**40**, **interaction.depth**=**8**, **shrinkage**=**0.09**, **n.minobsinnode**=**12**.

### Conclusiones

```{r}
results_grid_fine <- model_gbm_fine$results
results_grid_fine$Method <- "Fine Tune Grid"

# cominamos los resultados de los 3 modelos
combined_results <- bind_rows(results_auto, results_grid, results_grid_fine)

# Creamos el boxplot
ggplot(combined_results, aes(x = Method, y = Accuracy, fill = Method)) +
  geom_boxplot() +
  stat_summary(fun = mean, geom = "text", aes(label = round(after_stat(y), 3)), vjust = -0.5) +
  theme_minimal() +
  labs(
    title = "Comparación de Accuracy entre Métodos",
    x = "Método",
    y = "Accuracy",
    fill = "Método"
  )



```

Considerando los resultados finales, vemos que todos tienen un Accuracy muy parecido, aunque con el **ajuste fino** hemos subido la media un punto prácticamente. Además, al hacer **fine tuning** hemos **reducido considerablemente la variabilidad**. Esto último acompañado de que genera pocos valores atípicos, nos dice que este modelo daría es un **resutado más fiable**. 

#### Prediciendo nuevos valores y evaluación final del modelo

```{r}
# Extraemos la mejor fila de cada conjunto de resultados
best_auto <- results_auto[which.max(results_auto$Accuracy), ]
best_grid <- results_grid[which.max(results_grid$Accuracy), ]
best_grid_fine <- results_grid_fine[which.max(results_grid_fine$Accuracy), ]

# Agregamos una columna para identificar el método
best_auto$Method <- "Automático"
best_grid$Method <- "Grid Manual"
best_grid_fine$Method <- "Grid Fine Tuning"

# Combinar los mejores resultados y reordenamos las columnas
best_results <- rbind(best_auto, best_grid, best_grid_fine)
best_results <- select(best_results, Method, Accuracy, everything())

# Mostrar el data frame
print(best_results)

```

Viendo estos resultados podemos decir que, a priori, el **modelo con fine tuning** es el mejor de los 3. Podemos ver que obtiene mejor puntuación tanto en **accuracy** (aunque la diferencia es muy poca) como en **kappa**, lo que sugiere que es más consistente que los demás y además es más preciso. Además, comprobando tanto **accuracySD** como **kappaSD**, es el que menos desviación estándar tiene en ambos valores. Todo lo anterior refuerza nuestra hipótesis de que el modelo de **fine tuning** es el más estable, confiable y con mejor rendimiento de los 3.

Analizando las variables, podemos ver que el modelo de **fine tuning tiene un valor menor en shrinkage** (**0.09**), que le ayuda a tener más precisión. Su **interaction.depth** (**8**), el más alto junto con el grid manual, hace que pueda **detectar patrones más complejos**. También cuenta con el **n.minobsinnode** (**12**) más alto, **promoviendo la generalización**. Sin embargo, en cuanto a **n.trees** (**40**), es el más bajo. Esto, en este caso, **creemos que ayudaría a reducir la complejidad y mejorar la eficiencia**, pues el resto de parámetros confirman que es suficiente para no caer en **underfitting**.

## Calculando los valores sobre test de todos los modelos

Tras entrenar y evaluar los modelos sobre los conjuntos de datos, se realiza la comparación entre estos para decidir cual es el modelo **que mejor se adapta** al problema. Para ello, se usa el conjunto Test que se ha separado desde el inicio y se obtienen los resultados del modelo al predecir sobre este conjunto. Hay que tener en cuenta que, **para cada modelo**, **el conjunto de test usado habrá sufrido las mismas transformaciones que se aplicaron al conjunto de entrenamiento**.

### Test Binomial

Basarnos únicamente en el Accuracy dado tras predecir los modelos provocaría que estuviesemos incluyendo factores del azar, llevándonos a **conclusiones engañosas** puesto que no reflejaría el verdadero desempeño del modelo ante nuevos datos. Es por esta razón que se ha optado por realizar **pruebas estadísticas**, más concretamente, **tests binomiales**. Estas pruebas además de calcular el Accuracy del modelo, nos permite conocer un rango de valores (conocido como **intervalos de confianza**) en el que esta métrica se encuentra en la práctica.

Los intervalos de confianza nos permiten una **comparación sólida** entre los modelos **al evitar la aleatoriedad**. Además, podemos evaluar si los modelos han aprendido patrones generales relevantes de los datos si los comparamos con estrategias básicas como **predecir siempre hacia la clase mayoritaria**.

Para cada modelo se ha calculado mediante la siguiente función:
* **Accuracy** obtenido al predecir sobre el conjunto Test.

* **NIR**: proporción de la clase mayoritaría. 

* **lowerCI95**: Limite inferior del Intervalo de Confianza del 95%.

* **upperCI95**: Límite superior del Intervalo de Confianza del 95%.

* **pvalAcc>NIR**: Test binomial cuyo objetivo es determinar si el rendimiento del modelo entrenado es significativamente mejor que el rendimiento ofrecido con la predicción por clase mayoritaria.


```{r}
getAccModelList <- function(modelListWithData) {
  # modelListWithData es una lista donde cada elemento es una lista que incluye el modelo 
  # y el conjunto de datos correspondiente (preprocesado si es necesario).

  test_acc <- function(model, data, varSal, propMaxClass) {
    # Calculamos proporción de la clase mayoritaria (NIR)
    tbl <- table(data[[varSal]])
    propMaxClass <- max(tbl) / sum(tbl)

    # Predicción
    pred <- predict(model, data[, !names(data) %in% varSal, drop=FALSE])

    # Calculamos cuántos aciertos tenemos
    correctos <- sum(pred == data[[varSal]])
    n <- nrow(data)

    # Exactitud
    accuracy <- correctos / n

    # Test binomial para intervalo de confianza al 95%
    t1 <- binom.test(correctos, n, conf.level=0.95)
    lowerCI <- t1$conf.int[1]
    upperCI <- t1$conf.int[2]

    # Test binomial para comparar si Accuracy > NIR
    t2 <- binom.test(correctos, n, p = propMaxClass, alternative = "greater")
    pval_AccMayorNIR <- t2$p.value

    # Devolvemos vector con resultados
    c(lowerCI, accuracy, upperCI, pval_AccMayorNIR, propMaxClass)
  }

  # Aplicamos test_acc a cada modelo y su conjunto de datos correspondiente
  out <- lapply(modelListWithData, function(entry) {
    model <- entry$model
    data <- entry$data
    varSal <- entry$varSalida

    test_acc(model, data, varSal, propMaxClass = NULL)
  })

  # Convertimos la lista en data.frame
  out <- do.call("rbind", out)
  colnames(out) <- c("lowerCI95", "Accuracy", "upperCI95", "pvalAcc>NIR", "NIR")

  as.data.frame(out)
}

# Ejemplo de uso con conjuntos de datos preprocesados
listaModelosConDatos <- list(
  RPART = list(model = modelo_rpart, data = credit.Datos.Test, varSalida = "V16"),
  RF = list(model = modelo_rf, data = credit.Datos.Test, varSalida = "V16"),
  NNET = list(model = nnet_modelManual, data = credit_dummies_test_normalizado, varSalida = "V16"),
  GBM = list(model = model_gbm_fine, data = credit.Datos.Test, varSalida = "V16")
)

credit.results.test <- getAccModelList(listaModelosConDatos)
format(credit.results.test, digits=2, nsmall=3, scientific=F)


```


Una vez obtenido los resultados, podemos ver como los **intervalos de confianza** de todos los modelos **se superponen entre sí**. Esto nos quiere decir que **no hay una diferencia estadisticamente significativa** entre los modelos en cuanto a su **rendimiento**. El rendimiento de cada modelo está contenida en el intervalo de confianza, y **puede variar según la muestra**. Sin embargo, **Random Forest** presenta un intervalo un poco más ajustado y superior con respecto al resto, por lo tanto la variabilidad en cuanto al desempeño de **Random Forest** será menor a comparación de otros modelos, además de que puede alcanzar precisiones más altas que el resto de modelos, por lo que, **Random Forest es el modelo final que hemos elegido**.

Todos los modelos presentan unos intervalo de confianza que muestran un desempeño de los modelos aceptable. No difieren mucho del Accuracy obtenido en el entrenamiento y validación del modelo, lo que nos puede llevar a pensar que no sufren de **overfitting** ni **underfitting**. Esto tambien se domprueba con el campo **pvalAcc>NIR**. Todos los modelos tienen un valor significativamente bajo (siendo Random Forest el que presenta un menor valor), lo que muestra que estos modelos han aprendido patrones **significativos** y **generales** de los datos.


```{r}
plotTestRes <- function(x, ...) {
  plotCI(
    1:nrow(x),
    y = x$Accuracy,
    uiw = x$upperCI95 - x$Accuracy,
    ylab = "Accuracy",
    xlab = "Methods",
    xaxt = "n",
    ...
  )
  axis(1, at = 1:nrow(x), labels = row.names(x))
  grid (NULL, NULL, lty = 6, col = "cornsilk2")
  NIR <- x["NIR", 2]
  if (!is.na(NIR))
    abline(h = NIR, lty = 3)
}

# Graficar resultados
plotTestRes(credit.results.test)
```