main.Rmd

---
title: "Práctica de Aprendizaje Estadístico - Parte 1"
author: "Eduardo Valero Vilella y Cristian Calva Troya"
date: "2024-11-10"
output:
  html_document:
    toc: true
    toc_float: true
---
```{r}
options(warn = -1)
suppressPackageStartupMessages({
  library(plotly)
  library(dplyr)
  library(glmnet)
  library(corrplot)
  library(caret)
})
```

# 0. Carga los datos y elimina la variable TRAIN.

## 0.1 Configuración inicial. Carga de librerías.

**Nota: Nos aseguramos de tener el paquete `ElemStatLearn` instalado. Si no: `install.packages(lib_URL)`**

```{r}
lib_URL = 'https://cran.r-project.org/src/contrib/Archive/ElemStatLearn/ElemStatLearn_2015.6.26.2.tar.gz'
# install.packages(lib_URL)
# install.packages("plotly")

library(plotly)
library(ElemStatLearn)
library(dplyr)
```

Utilizamos `head()` para ver las primeras filas y hacernos una idea de las variables disponibles

```{r }
data("prostate", package = "ElemStatLearn")
head(prostate)
attach(prostate)
```

## 0.2 Eliminación de datos `TRAIN`.

La columna `train` se utiliza en estudios previos para indicar las filas de entrenamiento.

Dado que esta práctica no la requiere, eliminamos esta columna para simplificar el conjunto de datos.

```{r}
prostate = prostate %>% select(-train)
head(prostate)
```

# 1. Exploración de los datos

## 1.1 ¿Cuántas variables hay?

Tras eliminar la columna `train`, contamos con un total de `r length(prostate)` variables.

## 1.2 ¿De qué clase son?

Utilizamos la función `sapply` para aplicar una función, en este caso `class`, a todas las columnas.

```{r}
sapply(prostate, class)
```

Los tipos de las variables son del tipo `numeric` e `integer`, por lo que todos datos son números. Dado que son datos numéricos, podemos echar un vistazo a sus distribuciones.

```{r}
par(mfrow=c(3,3))
for (column in names(prostate)) {
  hist(prostate[[column]], main=paste("Histograma de", column), xlab = column)
}
```

Podemos observar que los valores de la variable `svi` se concentran en el 0 y en el 1, indicador de que esta variable podría ser un *booleano*. Veamos si estos son sus únicos valores.

```{r}
unique(prostate$svi)
```

0 y 1 son los únicos valores presentes, por lo que confirmaríamos que se trata de una variable *booleana* donde 1 indicaría que el cáncer ha invadido el vesículo seminal, y 0 indicaría lo contrario.

## 1.3 ¿Hay una variable que correspondiente al identificador de paciente? Si es así, elimínala.

No existe ninguna variable con el identificador del paciente.

## 1.4 ¿Hay valores nulos en alguno de los ficheros?

Averiguamos si existe algún valor nulo en todo el conjunto de datos. En caso de que lo haya, identificamos las columnas que los contienen.

```{r}
any(is.na(prostate))
```

Dado que no hay ningún valor nulo en nuestro set de datos, continuamos.

## 1.5 ¿Están estandarizadas las variables? En este punto del análisis, ¿es necesario normalizarlas?

Tras revisar los histogramas realizados en el subapartado 2 podemos concluir que, dados los rangos de las variables y su disposición, estos no se encuentran normalizados.

## 1.6 ¿Por qué crees que algunas variables están es escala logarítmica?

El objetivo de poner algunas de las variables en escala logarítmica es el de facilitar su posterior análisis por varios motivos:

-   Reducir la asimetría en casos en los que la distribución está muy sobrecargada en uno de sus extremos donde algunos valores excepcionales distorsionarían las estadísticas de la variable, haciendo que sus distribuciones se asemejen más a las de una distribución normal (como se observa en `lcavol`, `lweight` y `lpsa`).
-   Hacer que algunas variables se comporten de forma lineal cuando en escala natural no lo harían.

Además, cabe añadir que en los campos de la biología y la química es común usar escalas logarítmicas cuando estas hacen más fácil la interpretación de los datos.

# 2. Análisis de variables categóricas

Para convertir las variables `svi` y `gleason` a categóricas usamos la función `factor`.

```{r}
prostate$svi <- factor(prostate$svi)
prostate$gleason <- factor(prostate$gleason)
```

A continuación, visualizamos las distribuciones usando gráficos de barras:

```{r}
par(mfrow = c(1, 2))  # Ajustar disposición gráfica
barplot(table(prostate$svi), main = "Distribución de SVI", col = "lightblue")
barplot(table(prostate$gleason), main = "Distribución de Gleason", col = "lightgreen")
```

Se puede observar como para la mayoría de pacientes se indica que el cáncer no ha invadido la vesícula seminal, a la par que la mayoría tiene un índice de Gleason de 6 o 7. Esto puede indicar cierta relación, ya que el índice de Gleason mide la agresividad del cáncer, por lo que índices más bajos indican cánceres menos agresivos. Podríamos aventurarnos a pensar que la mayoría de cánceres agresivos no llegarían a infectar la vesícula seminal.

En el caso de la variable **Edad** queremos hacer una agrupación cada 5 años, por lo que usaremos la función `cut` para generar esas categorías:

```{r}
max.age = max(prostate$age)
min.age = min(prostate$age)
# Añadimos +5 al maximo para que el rango tome tambien el ultimo grupo
prostate$age <- cut(prostate$age, breaks = seq(from=min.age, to=max.age+5, by=5), right=FALSE)
barplot(table(prostate$age), main = "Distribución de Edad", col = "salmon", width = 3)
```

En el caso de la edad nos encontramos una distribución que podría recordar a la normal, centrada en los 65 años, encontrándose el grueso de la distribución entre los 61 y 70 años. Cabe añadir que encontramos casos desde los `r min.age` hasta los `r max.age` años.

# 3. Análisis de frecuencias

## 3.1 ¿Qué porcentaje de pacientes con la puntuación de Gleason igual a 7, presenta índice igual svi igual a 0?

```{r}
# Filtramos por gleason = 7
gleason.7 <- subset(prostate, gleason == 7)
gleason.7.props <- prop.table(table(gleason.7$svi))*100
```

Por lo tanto, un `r round(gleason.7.props["0"])`% de los pacientes con índice de Gleason igual a 7 presenta un `svi` de 0. Podemos representar esta distribución mediante un gráfico de sectores:

```{r}
fig <- plot_ly(
  labels = names(gleason.7.props),
  values = as.numeric(gleason.7.props),
  type="pie",
  textinfo = "label+percent",
  hoverinfo = "none",
  insidetextorientation = "radial",
  height = 400
) %>%
  layout(
    title = "Distribución del índice de Gleason para pacientes con svi = 0",
    showlegend = TRUE
  )
fig
```

## 3.2 ¿Qué porcentaje de pacientes con índice svi igual a 0 tiene la puntuación de Gleason igual a 7?

```{r}
# Filtramos por gleason = 7
svi.0 <- subset(prostate, svi == 0)
svi.0.props <- prop.table(table(svi.0$gleason))*100
```

Por lo tanto, un `r round(svi.0.props["7"])`% de los pacientes con un `svi` de 0 presenta un índice de Gleason igual a 7. Podemos visualizar la distribución del índice de Gleason para pacientes con un `svi`$=0$ mediante un gráfico de sectores:

```{r}
fig <- plot_ly(
  labels = names(svi.0.props),
  values = as.numeric(svi.0.props),
  type="pie",
  textinfo = "label+percent",
  hoverinfo = "none",
  insidetextorientation = "radial" ,
  height = 400
) %>%
  layout(
    title = "Distribución del índice de Gleason para pacientes con svi = 0",
    showlegend = TRUE
  )
fig
```

## 3.3 Estas dos variables, ¿son independientes?

Vistos los resultados de estos dos últimos subapartados, podemos pensar que estas variables están relacionadas. Además, ya en el apartado 2 se vio como también podría existir una posible relación entre las variables y se le trató de dar una explicación lógica.

Para probar si realmente existe esta relación entre variables, exponemos los datos a un test de independencia *chi cuadrado*, cuya implementación en R es simple y rápida. En este test queremos probar que las variables están relacionadas, lo cual es nuestra hipótesis $H_1$, por lo que la hipótesis nula $H_0$ será que estas son independientes:

```{r}
contingencia.svi_gleason <- table(prostate$gleason, prostate$svi)
chi_test <- chisq.test(contingencia.svi_gleason)
chi_test
```

Al realizar el test, obtenemos un $p$-valor de 0.0012. Dado que los niveles de significación más comunes son $0.05$ y $0.01$, y nuestro $p$-valor es notablemente menor, podemos afirmar con seguridad que existe una relación entre estas variables.

# 4. Regresión lineal simple.

Análisis de la dependencia lineal de la variable `lpsa` con la variable `lcavol`.

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

Ajustamos el modelo y vemos su `summary()`:

```{r}
modelo_lineal <- lm(lpsa ~ lcavol, data = prostate)

summary(modelo_lineal)
```

## 4.1 Interpretación del Modelo Lineal calculado

Ecuación del Modelo:

La ecuación del modelo de regresión lineal es:

$$
\text{lpsa} = 1.50730 + 0.71932 \cdot \text{lcavol}
$$

Intercepto ($\beta_0$):

-   El intercepto de **1.50730** indica el valor estimado de `lpsa` cuando `lcavol` es igual a 0.

Pendiente ($\beta_1$):

-   El coeficiente para `lcavol` es **0.71932**. Esto sugiere que, por cada unidad de incremento en `lcavol`, se espera que `lpsa` aumente en aproximadamente **0.71932** unidades.

## 4.2 Plot de los datos junto a la recta de regresión.

Crearemos una gráfico de dispersión de `lcavol` vs `lpsa`.

```{r}
ggplot(prostate, 
       aes(x = lcavol, y = lpsa)) +
  
 # Puntos de los datos
geom_point(color = "blue", 
           alpha = 0.4) +

# Recta de regresión
geom_smooth(method = "lm", 
            color = "red") + 
  
labs(title = "Relación entre lcavol y lpsa",
     x = "Log de Volumen de Cáncer (lcavol)",
     y = "Log de Antígeno Prostático Específico (lpsa)") +
  
theme_minimal()
```

## 4.3 Intervalos de confianza para los coeficientes del modelo con una confianza de 0.95

Los intervalos de confianza al 95% para los coeficientes del modelo son importantes para entender la precisión de las estimaciones. A continuación se presentan los intervalos calculados:

```{r}
confint(modelo_lineal, level = 0.95)
```

Con esto podemos estimar que con ausencia de un tumor medido (`lcavol`), el `psa` se encontrará entre valores de 1.26 y 1.74 ng/m, en su medida logarítmica.

Ademas, podemos estimar que por cada cm (10mm) de tumor aumentado, el PSA aumenta en unidades comprendidos entre 5.8 y 8.5 ng/m.

En el análisis de regresión lineal, se ha observado que tanto el coeficiente del intercepto como el coeficiente de `lcavol` presentan valores $p$ muy pequeños $(< 2e-16)$. Esto indica que podemos rechazar la hipótesis nula $(H_0)$ de que los "verdaderos" coeficientes son nulos. Específicamente, esto significa que existe una relación estadísticamente significativa entre el volumen del cáncer (`lcavol`) y el nivel de antígeno específico de la próstata (`lpsa`).


## 4.4 RSE y la eficacia del modelo.

El **Error Estándar de los Residuos (RSE)** es una medida de la cantidad de variabilidad de los residuos (la diferencia entre los valores observados y los valores predichos por el modelo) en un modelo de regresión lineal. Se utiliza para evaluar la precisión de las predicciones del modelo. Un RSE más bajo indica que el modelo tiene un mejor ajuste a los datos.

Matemáticamente, el RSE se define como:

$$
RSE = \sqrt{\frac{1}{n-2} \sum (y_i - \hat{y}_i)^2}
$$

donde:

-   $y_i$ son los valores observados,

-   $\hat{y}_i$ son los valores predichos por el modelo,

-   $n$ es el número de observaciones.

Si echamos un vistazo de nuevo a nuestro resultado:

```{r}
summary(modelo_lineal)
```

### RSE

Un RSE de 0.7875 indica que, en promedio, las predicciones del modelo se desvían de los valores observados de la variable de respuesta (`lpsa`) en aproximadamente 0.7875 unidades en la escala logarítmica.

Indicando que por cada cm(10mm) de aumento en el `lcavol` tenemos un desvía de 7.8 ng/m.

### R-cuadrado (R²)

El coeficiente de determinación $R^2$ es una métrica que indica la proporción de la varianza en la variable dependiente `lpsa` que es predecible a partir de la variable independiente `lcavol`. En este modelo, el $R^2$ es:

$R^2 = 0.5394$

-   **Interpretación**: Aproximadamente el 53.94% de la variabilidad en los niveles de PSA puede ser explicada por el volumen de cáncer. Esto sugiere una relación moderada entre las variables.

### R-cuadrado Ajustado (Adjusted R²)

El $R^2$ ajustado, que en este caso es:

$R^2$ Ajustado = 0.5346

-   **Interpretación**: El $R^2$ ajustado penaliza el $R^2$ por el número de predictores en el modelo, lo que lo hace más adecuado cuando se comparan modelos con diferentes números de variables. Este valor indica que el modelo tiene un ajuste similar al $R^2$, sugiriendo que el volumen de cáncer sigue siendo un predictor significativo.

## 4. 6 Interpretación del Modelo Lineal (`psa` y `cavol`)

### Ecuación del Modelo

El modelo de regresión lineal simple que hemos calculado es:

$$ 
\text{lpsa} = 1.50730 + 0.71932 \times \text{lcavol} 
$$

### Relación entre PSA y cavol

1.  **Modelo Original**: La ecuación del modelo de regresión lineal que obtuvimos es:

  $$
    \text{lpsa} = 1.50730 + 0.71932 \times \text{lcavol}
  $$

  $$
    ln(\text{psa}) = 1.50730 + 0.71932 \times ln(\text{cavol})
  $$

2.  **Despejar PSA**: Como $\text{lpsa}$ es el logaritmo natural del PSA, puedes deshacerte del logaritmo aplicando la función exponencial:

  $$
    \text{PSA} = e^{ln(\text{psa})} = e^{(1.50730 + 0.71932 \times ln(\text{cavol}))}
  $$

  $$
    \text{PSA} = e^{1.50730 } + e^{(0.71932 \times ln(\text{cavol}))}
  $$

  $$
    \text{PSA} = e^{1.50730 } + (e^{ln(\text{cavol})})^{0.71932}
  $$

3.  **Transformación a `cavol`**: Para expresar esto en términos de `cavol` (el volumen de cáncer en su forma original), tienes que revertir el logaritmo:

  $$
    \text{cavol} = e^{ln(\text{cavol})}
  $$

Ahora sustituimos `lcavol` por su expresión en términos de `cavol`:

  $$
    \text{PSA} = e^{1.50730 } +\text{cavol}^{0.71932} 
  $$

4.  **Relación Final**: La relación entre `PSA` y `cavol` queda como:

  $$
    \text{PSA} = C \cdot \text{cavol}^{b}
  $$ 
    Donde $C = e^{1.50730}$ y $b = 0.71932$. Esto significa que `PSA` es proporcional a `cavol` elevado a un exponente.

### Interpretación de la Relación

-   **C**: Este es un coeficiente constante que representa el nivel de `PSA` cuando `cavol` es igual a 1 (o la unidad de referencia).
-   **b (0.71932)**: Indica que por cada aumento en el volumen de cáncer (`cavol`), el `PSA` aumenta en un porcentaje específico que depende del valor de `cavol`.

# 5. Regresión lineal múltiple.

## 5.1 Observación de correlaciones entre variables.

```{r}
library(corrplot)
prostate$svi = as.numeric(prostate$svi)

prostate2 = prostate

prostate2$age_num <- as.numeric(prostate2$age)
prostate2$pgg45_num <- as.numeric(prostate2$pgg45)
prostate2$gleason_num <- as.numeric(prostate2$gleason)

cor_data = prostate2 %>% select(-age, -pgg45, -gleason)
cor_matrix = cor(cor_data)

corrplot(cor_matrix, 
         diag = FALSE,
         method = "circle", 
         type = "upper", 
         tl.col = "black", 
         tl.srt = 25,
         addCoef.col = "black")
```

La matriz de correlación muestra las correlaciones de Pearson entre las variables del conjunto de datos, con valores que van de -1 a 1.

### Análisis de la Matriz de Correlación

Se ha realizado un análisis de correlación entre las variables del conjunto de datos. A continuación se presentan los resultados y las interpretaciones:

1.  **(`lcavol`, `lpsa`)**
    -   **Máxima correlación mutua**: (0.734).
    -   **Interpretación**: La fuerte correlación entre el volumen del cáncer en escala logarítmica (`lcavol`) y el nivel de antígeno específico de la próstata (PSA) sugiere que a medida que aumenta el volumen del cáncer, también tienden a aumentar los niveles de PSA.
2.  **(`lbph`, `lweight`)**:
    -   **Máxima correlación mutua**: (0.442).
    -   **Interpretación**: La correlación entre el peso en escala logarítmica (`lweight`) y la hipertrofia benigna de próstata (`lbph`) indica que hay una relación positiva entre ambos. Esto sugiere que los hombres con mayor peso pueden tener mayor riesgo de desarrollar condiciones benignas en la próstata.
3.  **`lcp`**:
    -   **Máxima correlación**: `lcavol` (0.675).
    -   **Interpretación**: La correlación entre el nivel de extensión del cáncer dentro de la cápsula (`lcp`) y el volumen del cáncer (`lcavol`) tienen trivial relación.
4.  **`svi`**:
    -   **Máxima correlación**: `lcp` (0.68).
    -   **Interpretación**: La correlación entre el nivel de invasión seminal (`svi`) y el nivel de extensión del cáncer dentro de la cápsula (`lcp`) es positiva y moderada (0.68). Esto indica que a medida que aumenta la invasión a las vesículas seminales, también podría aumentar la extensión del cáncer dentro de la cápsula.
5.  **`age_num`**:
    -   **Máxima correlación:** `lbph` (0.364).
    -   **Interpretación:** La correlación moderada entre la edad en formato numérico (`age_num`) y la hipertrofia prostática benigna (`lbph`) sugiere que a mayor edad, la probabilidad de presentar hipertrofia benigna aumenta.
6.  **(`pgg45_num`, `gleason_num`)**:
    -   **Máxima correlación mutua:**  (0.752).
    -   **Interpretación:** La alta correlación entre el porcentaje de patrón Gleason 4 o 5 (`pgg45_num`) y la puntuación total de Gleason (`gleason_num`) refuerza que un mayor porcentaje de patrones más agresivos se asocia con una puntuación de Gleason más alta.

## 5.2 Modelo de regresión lineal para predecir `lpsa`.

```{r}
modelo_lpsa = lm(lpsa ~ lcavol + lweight + lbph + lcp + svi, 
                  data = prostate)

summary(modelo_lpsa)
```

Siendo la siguiente su función:

$$
\text{lpsa} = -1.12471 + 0.54790 \cdot \text{lcavol} + 0.53036 \cdot \text{lweight} + 0.07999 \cdot \text{lbph} - 0.03638 \cdot \text{lcp} + 0.75975 \cdot \text{svi} + \epsilon
$$

### Coeficientes:

-   `lcavol`: 0.54790 (muy significativo, `p = 8.56e-09`). Esto indica que, manteniendo las otras variables constantes, un aumento de una unidad en el volumen del cáncer (`lcavol`) se asocia con un aumento de 0.54790 unidades en los niveles de PSA (`lpsa`).

-   `lweight`: 0.53036 (significativo, `p = 0.00867`). Similar a `lcavol`, un aumento de una unidad en el peso del paciente (`lweight`) se asocia con un aumento en los niveles de PSA, sugiriendo que el peso del paciente podría influir en el diagnóstico.

-   `lbph`: 0.07999 (no significativo, `p = 0.15971`). La hipertrofia benigna de próstata (`lbph`) no tiene un efecto significativo en los niveles de PSA, ya que su valor de p es mayor que 0.05.

-   `lcp`: -0.03638 (no significativo, `p = 0.65391`). El nivel de extensión del cáncer dentro de la cápsula (`lcp`) no tiene un efecto significativo sobre los niveles de PSA.

-   `svi`: 0.75975 (muy significativo, `p = 0.00221`). La invasión seminal (`svi`) tiene una asociación significativa con los niveles de PSA, indicando que a medida que la invasión seminal aumenta, también aumentan los niveles de PSA.

### Rendimiento

$R^{2}$ ajustado: 0.6248. El modelo explica aproximadamente el 62.48% de la variación en los niveles de PSA. Esto sugiere una mejora en la capacidad predictiva con la adición de la variable svi, aunque aún hay un 37.52% de la variación que no se explica por las variables incluidas en el modelo.

**RSE (Error estándar de los residuos):** 0.7071. Este valor indica el tamaño promedio de los errores en las predicciones del modelo, siendo relativamente bajo, lo que sugiere que las predicciones son bastante precisas.

**F-estadístico** 32.97 con p-valor \< 2.2e-16. Este resultado indica que el modelo es altamente significativo en su conjunto, lo que significa que al menos una de las variables predictoras tiene un impacto significativo en los niveles de PSA.

### Conclusiones

El modelo de regresión muestra que `lcavol`, `lweight`, y `svi` son los factores más relevantes para predecir los niveles de PSA, con relaciones significativas y positivas. En cambio, `lbph` y `lcp` no tienen un impacto significativo en los niveles de PSA. El $R^{2}$ ajustado de 0.6248 indica que el modelo explica el 62.48% de la variabilidad en PSA, y el RSE bajo sugiere predicciones precisas. El modelo es útil, pero podría mejorarse considerando otro conjunto de variables.

# 6. Modelo de Ridge y Lasso

Para realizar las predicciones usando los modelos de **Ridge** y **Lasso** utilizaremos la librería `glmnet`. Primero preparamos nuestro conjunto de datos lo dividimos entre el conjunto de **training** y el de **testeo**. Para crear nuestro training set utilizaremos el 70% de las entradas de nuestro conjunto de datos.

```{r}
library(glmnet)

# Preparamos un dataset sin las variables que no usamos
datos_usados <- subset(prostate, select = -c(age, pgg45, gleason))

# Dividimos el dataset entre predictores y el valor a predecir
X <- as.matrix(datos_usados[, -which(names(datos_usados) == "lpsa")])
y <- datos_usados$lpsa

# Creamos los conjuntos de train y test, especificando la seed para poder reproducir la ejecución
# con los mismos sets
set.seed(123) 
train_idx <- sample(1:nrow(X), size = 0.7 * nrow(X))
X.train <- X[train_idx,]
y.train <- y[train_idx]
X.test <- X[-train_idx,]
y.test <- y[-train_idx]
```

## 6.1 Modelo de Ridge

Para utilizar el modelo de Ridge con `glmnet` debemos poner el parámetro $\alpha = 0$. Utilizaremos un *grid* con 200 valores de $\lambda$ diferentes, desde $10^{4}$ hasta $10^{-4}$. Estos valores se han determinado después de realizar la gráfica varias veces para asegurarnos de que las regiones que en ella se observan tienen una representación similar. A partir del modelo, podemos obtener el conjunto de coeficientes y los lambdas correspondientes a cada uno de los 200 modelos obtenidos.

```{r}
lambdas = 10^seq(4,-4, length = 200)

ridge.model <- glmnet(X.train, y.train, alpha = 0, lambda = lambdas)
# Obtenemos los coeficientes y los lambdas
ridge.coefs <- as.matrix(ridge.model$beta)
ridge.lambda <- ridge.model$lambda

matplot(log(ridge.lambda), t(ridge.coefs), type = "l", lty = 1, col = rainbow(nrow(ridge.coefs)),
        xlab = "log(λ)", ylab = "Coeficientes", main = "Ridge: Coeficientes vs log(λ)")
legend("topright", legend = colnames(X), col = rainbow(nrow(ridge.coefs)), lty = 1, cex = 0.6)
```

Vemos claramente 3 regiones en la gráfica. Una primera región a la izquierda, hasta $\log\lambda \approx -4$, donde los coeficientes son constantes. Una segunda región intermedia, aproximadamente entre $-4$ y $4$ donde estos valores varían, disminuyendo todos los coeficientes a excepción de `lcp`, el cual comienza con un valor negativo y crece hasta tener uno positivo. Finaliza con una tercera región, donde los coeficientes convergen a un valor cercano a 0. 

Obtenemos el valor óptimo estimado de $\lambda$ usando validación cruzada, con la función `cv.glmnet` con el parámetro $\alpha=0$:

```{r}
# La seleccion hecha en el CV es tambien aleatoria, asi que debebemos
# especificar la seed para asegurarnos de que podemos reprocir la prueba
set.seed(123) 
# Obtenemos el mejor valor de lambda usando cross-validation
cv.ridge.model <- cv.glmnet(X.train, y.train, alpha=0, lambda=lambdas)
lambda.min.ridge <- cv.ridge.model$lambda.min
```

- Valor óptimo: $\lambda_\text{Ridge} = `r round(lambda.min.ridge,4)`$.
- Valor óptimo: $\log\lambda_\text{Ridge} = `r round(log(lambda.min.ridge),4)`$.

El valor óptimo del $\log\lambda_\text{Ridge}$ nos permite ubicar el valor del $\lambda$ óptimo en la gráfica anterior, ubicándola en la región intermedia descrita con anterioridad.

A continuación presentamos la gráfica del $MSE$ frente a $\log\lambda$ correspondiente al modelo de Ridge.

```{r}
plot(cv.ridge.model)
```

Vamos ahora a hacer la predicción de `lpsa` usando el conjunto de testeo.

```{r}
pred.ridge <- predict(cv.ridge.model, s = lambda.min.ridge, newx = X.test)
```

## 6.2 Modelo de Lasso

Para utilizar el modelo de Lasso con `glmnet` debemos poner el parámetro $\alpha = 1$. Utilizaremos el mismo *grid* de \lambda usado para el modelo de *Ridge*, pero con valores de $\lambda$ entre $10^1$ y $10^{-4}$, valores que obtenemos de igual manera tras reproducir la gráfica varias veces. Después procedemos de forma similar.

```{r}
lambdas = 10^seq(1,-4, length = 200)
lasso.model <- glmnet(X.train, y.train, alpha = 1, lambda = lambdas)
# Obtenemos los coeficientes y los lambdas
lasso.coefs <- as.matrix(lasso.model$beta)
lasso.lambda <- lasso.model$lambda

matplot(log(lasso.lambda), t(lasso.coefs), type = "l", lty = 1, col = rainbow(nrow(lasso.coefs)),
        xlab = "log(λ)", ylab = "Coeficientes", main = "Lasso: Coeficientes vs log(λ)")
legend("topright", legend = colnames(X), col = rainbow(nrow(lasso.coefs)), lty = 1, cex = 0.6)
```

Podemos observar un gráfico similar al obtenido para *Ridge*, con la diferencia de que *Lasso* sí elimina predictores, es decir hace sus coeficientes 0. Contando esta diferencia, vemos también una representación análoga de las 3 regiones que vimos en el gráfico del modelo de *Ridge*, con una primera región donde el coeficiente de `lcp` es negativo, hasta que se elimina, una región intermedia donde se regulan el resto de parámetros, y una final donde se eliminan todos los coeficientes (en *Ridge* también tendían a 0). Es por esto que podemos suponer que, también en este caso, el valor óptimo de $\lambda$ se encontrará en la región intermedia.

También podemos visualizar el `plot` por defecto del modelo de *Lasso*, el cual nos muestra la evolución de los coeficientes con el `L1 Norm`, que representa el error que se intenta minimizar al utilizar este modelo:

$$
L1_\text{Norm} = RSS + \lambda \sum ^p_{j=1}|\beta_j|
$$

```{r}
plot(lasso.model, col = rainbow(nrow(lasso.coefs)))
legend("topright", legend = colnames(X), col = rainbow(nrow(lasso.coefs),start=0), lty = 1, cex = 0.6)
```

Obtenemos el valor óptimo estimado de $\lambda$ usando validación cruzada, con la función `cv.glmnet` con el parámetro $\alpha=1$:

```{r}
# La seleccion hecha en el CV es tambien aleatoria, asi que debebemos especificar la 
# seed para asegurarnos de que podemos reprocir la prueba
set.seed(123)
# Obtenemos el mejor valor de lambda usando cross-validation
cv.lasso.model <- cv.glmnet(X.train, y.train, alpha=1, lambda=lambdas)
lambda.min.lasso <- cv.lasso.model$lambda.min
log(lambda.min.lasso)
```

- Valor óptimo: $\lambda_\text{Ridge} = `r round(lambda.min.lasso,4)`$.
- Valor óptimo: $\log\lambda_\text{Ridge} = `r round(log(lambda.min.lasso),4)`$.

A continuación presentamos la gráfica del $MSE$ frente a $\log\lambda$ correspondiente al modelo de Lasso.

```{r}
plot(cv.lasso.model)
```

Vamos ahora a hacer la predicción de `lpsa` usando el conjunto de testeo.

```{r}
pred.lasso <- predict(cv.lasso.model, s = lambda.min.lasso, newx = X.test)
```

## 6.3 Comparación con modelo lineal múltiple

1. **Coeficientes**: Comparamos los coeficientes obtenidos en cada modelo.

Creamos un `data.frame` que contenga los coeficientes para cada modelo. Esta tabla nos sirve tanto para poder comparar cuantitativamente los valores de cada coeficiente, como para preparar los datos para presentarlos en un gráfico usando `ggplot`.
```{r}
coef.multi <- coef(modelo_lpsa)[-1]
coef.ridge <- coef(cv.ridge.model, s=lambda.min.ridge)[-1]
coef.lasso <- coef(cv.lasso.model, s=lambda.min.lasso)[-1]

coef.data <- data.frame(
  variable = names(coef.multi),
  multiple = as.numeric(coef.multi),
  ridge = as.numeric(coef.ridge),
  lasso = as.numeric(coef.lasso)
)

coef.data
```

Cambiamos la tabla de datos a formato `wide` para poder mapearlo correctamente usando ggplot

```{r}
coef.data.long <- reshape2::melt(coef.data, id.vars = "variable", 
                                 variable.name = "Modelo", value.name = "Coeficiente")

ggplot(coef.data.long, aes(x = variable, y = Coeficiente, fill = Modelo)) +
  geom_bar(stat = "identity", position = "dodge") + 
  labs(title = "Comparación de Coeficientes entre Modelos de Regresión",
    x = "Predictor", y = "Peso del Coeficiente") +
  theme_minimal()
```

- **`lbph`**: Lo primero que salta a la vista es que el modelo **Lasso** ha eliminado este predictor a diferencia del resto, ya que el modelo **Lasso** es el único que realiza selección de variables. De forma similar, el resto de modelos le otorgan un peso muy bajo.

- **`lcavol`**: Los tres modelos coinciden en otorgar un gran peso a este predictor, destacando el modelo lineal. Lasso y Ridge otorgan un valor similar entre ellos, siendo el peso en el modelo de Lasso mayor.

- **`lcp`**: Este caso es especialmente llamativo, y es que mientras los modelos Lasso y Ridge coinciden en otorgarle un peso alto (es el coeficiente de mayor valor en ambos), en el modelo de regresión múltiple este valor no es solo mucho menor, sino que además es negativo. 

- **`lweight`**: En este predictor encontramos el mayor consenso entre los tres modelos, otorgándole un gran peso los tres modelos.

- **`svi`**: Este es otro caso de gran contraste entre los modelos, y es que mientras es el valor con mayor peso en el modelo múltiple, con una gran distancia al segundo, el modelo de Ridge le otorga un valor mínimo y el modelo de Lasso directamente lo elimina.

El comportamiento tan diferente de los coeficientes `lcp` y `svi` podría indicar la existencia de una correlación de importancia notable entre dichas variables.

2. **Rendimiento**: 

Vamos a calcular la predicción del modelo lineal (el resto ya han sido calculados), y a continuación calcularemos y compararemos el **MSE** (*error cuadrático medio*) que resulta de cada modelo.

```{r}
# Funciones para calcular las métricas
mse <- function(y,y2) {
  se = (y-y2)^2
  mean(se)
}

rsq <- function(y,y2) {
  rss = sum((y-y2)^2)
  tss = sum((mean(y)-y)^2)
  1 - rss / tss
}

rsq.adj <- function(y,y2,p) {
  n = length(y)
  1 - ((1-rsq(y,y2))*(n-1))/(n-p-1)
}
```


```{r}
pred.multiple <- predict(modelo_lpsa, newdata = as.data.frame(X.test))

resultados <- data.frame(
  Metric = c("MSE", "R2", "R2adj"),
  Multiple = c(
    mse(y.test, pred.multiple), 
    rsq(y.test, pred.multiple), 
    rsq.adj(y.test, pred.multiple, 5)
  ),
  Ridge = c(
    mse(y.test, pred.ridge), 
    rsq(y.test, pred.ridge), 
    rsq.adj(y.test, pred.ridge, 5)
  ), 
  Lasso = c(
    mse(y.test, pred.lasso),
    rsq(y.test, pred.lasso),
    rsq.adj(y.test, pred.lasso, 4) # p = 4 aquí ya que, como recordamos, lasso elimina lbph
  )
)

resultados
```

Replicamos lo hecho para los coeficientes para graficar los resultados obtenidos:
```{r}
resultados.long <- reshape2::melt(resultados, id.vars = "Metric", variable.name = "Modelo", value.name = "Valor")

ggplot(resultados.long, aes(Metric, Valor, fill = Modelo)) +
  geom_bar(stat= "identity", position="dodge") +
  labs(title = "Comparación de Métricas de Rendimiento entre Modelos",
    x = "Métrica", y = "Valor de la Métrica") +
  theme_minimal()
```

## 6.4 Conclusiones

Como podemos observar, el modelo que mejor ha funcionado en nuestro caso es el de la regresión múltiple, teniendo un resultado mejor que el de los otros modelos tras medir su rendimiento usando **MSE**, $R^2$ y $R^2_\text{adj}$ en el conjunto de testeo.

Así, podemos concluir que los modelos que penalizan un número mayor de predictores, en este caso, funciona peor que el modelo de regresión múltiple que no cuenta con esa penalización. Además, tampoco representan una mejora en el rendimiento o en el espacio de almacenamiento, dado que contamos con un número muy limitado de predictores y de observaciones, por lo que no se recomendaría su uso en situaciones similares.

Para responder a si los modelos lineales podrían ser utilizados o no, revisamos los valores de las métricas obtenidas para el modelo de regresión múltiple:

- $R^2 = `r resultados$Multiple[2]`$
- $R^2_\text{adj} = `r resultados$Multiple[3]`$

Como vemos, nuestro modelo funciona relativamente bien con el conjunto de testeo, por lo que podemos concluir que este modelo podría llegar a ser utilizado para realizar predicciones.

# 7. LDA

```{r}
suppressPackageStartupMessages({
  library(MASS)
})
```

Ajustar el modelo LDA utilizando `svi` como variable de respuesta y `lcavol`, `lcp`, `lpsa` como variables predictoras.

```{r}
prostate$svi <- as.factor(prostate$svi)
lda_model <- lda(svi ~ lcavol + lcp + lpsa, data = prostate)
```

## 7.1 Visualizar los Coeficientes del Modelo

-   El modelo muestra las **probabilidades** de cada grupo (0 y 1). Estas indican la proporción de datos que pertenece a cada clase.

-   Las **medias de grupo** presentan los valores medios de las variables predictoras (`lcavol`, `lcp`, `lpsa`) para cada clase del objetivo `svi`.

-   Los **coeficientes de los discriminantes lineales** reflejan la importancia de cada variable en la construcción del discriminante lineal.

```{r}
lda_model
```

## 7.2 Predicciones con el Modelo LDA

-   `class`: contiene las clases predichas por el modelo.

-   `posterior`: contiene las probabilidades posteriores de cada clase para cada observación.

-   `x`: valores de los discriminantes lineales para cada observación.

```{r}
lda_model.pred <- predict(lda_model)
names(lda_model.pred)
```

## 7.3 Observar las Clases Predichas

Las primeras tres observaciones fueron clasificadas en la clase 0, lo que indica que el modelo ha predicho la ausencia de `svi` para estos casos.

```{r}
lda_model.pred$class[1:3]
```

## 7.4 Matriz de Confusión

-   La matriz de confusión muestra que el modelo clasificó correctamente 69 observaciones del grupo 0 y 17 del grupo 1.

-   Hay 7 casos del grupo 1 que fueron clasificados incorrectamente como 0, y 4 casos del grupo 0 clasificados incorrectamente como 1.

Esto podría estar relacionado con un **desbalance en las clases**, ya que el grupo 0 tiene una representación mucho mayor (aproximadamente el 78%) que el grupo 1, lo que puede hacer que el modelo esté más inclinado a predecir el grupo mayoritario.

Este desbalance puede afectar la capacidad del modelo para detectar correctamente los casos del grupo minoritario, lo que es crucial en escenarios donde los casos de la clase menos frecuente son más importantes, como en la detección de enfermedades.

```{r}
library(plotly)

conf_matrix <- table(lda_model.pred$class, svi)

conf_matrix_df <- as.data.frame(as.table(conf_matrix))

heatmap <- plot_ly(
  data = conf_matrix_df,
  x = ~Var1,
  y = ~svi,
  z = ~Freq,
  type = "heatmap",
  colors = c("white", "blue"),
  showscale = FALSE
) %>%
  layout(
    title = "Predicciones de SVI",
    xaxis = list(title = "Predicción"),
    yaxis = list(title = "Clase Real")
  )

heatmap %>%
  add_annotations(
    x = conf_matrix_df$Var1,
    y = conf_matrix_df$svi,
    text = conf_matrix_df$Freq, 
    showarrow = FALSE,
    font = list(size = 12, color = "black")
  )

```

## 7.5 Probabilidades Posteriores.

-   Las probabilidades posteriores indican cuán probable es que cada observación pertenezca a cada una de las clases (0 o 1).

-   En este caso, para las primeras tres observaciones, la probabilidad de pertenecer a la clase 0 es extremadamente alta (cerca de 1).

```{r}
lda_model.pred$posterior[1:3, ]
```

## 7.6 Valores de los Discriminantes Lineales

-   Los valores de los discriminantes lineales representan las puntuaciones de cada observación en la función discriminante.

-   Estos valores son útiles para entender cómo el modelo discrimina entre los grupos basándose en las variables predictoras.

```{r}
lda_model.pred$x[1:3]
```

## 7.7 Conclusiones

El modelo LDA ha mostrado un rendimiento desigual debido a un posible desbalanceo de clases, donde el grupo 0 es significativamente más grande que el grupo 1. Esto ha llevado a que el modelo clasifique correctamente 69 observaciones del grupo 0, pero solo 17 del grupo 1. Aunque el número de falsos positivos es relativamente bajo (4), los falsos negativos (7) indican que el modelo no está detectando adecuadamente todas las observaciones del grupo 1. Este comportamiento es común en escenarios con clases desbalanceadas.

Aunque el modelo no es completamente ineficaz, su capacidad para identificar el grupo minoritario es limitada. Se podrían explorar técnicas como el ajuste el balanceo de clases (por ejemplo, sobremuestreo o submuestreo) para mejorar el rendimiento en la clasificación del grupo 1 y obtener una predicción más equilibrada.


# 8. Regresión Logística

## 8.1 `svi` en función de `lcavol`, `lcp` y `lpsa`

Veamos primero las distribuciones del  `svi` frente a los predictores propuestos:

```{r}
boxplot1 <- plot_ly(data = prostate, x = ~factor(svi), y = ~lcavol, 
                    type = "box", name="lcavol", width = 900) %>%
  layout(xaxis = list(title = "svi"), yaxis = list(title = "lcavol"))

boxplot2 <- plot_ly(data = prostate, x = ~factor(svi), y = ~lcp, 
                    type = "box", name="lcp", width = 900) %>%
  layout(xaxis = list(title = "svi"), yaxis = list(title = "lcp"))

boxplot3 <- plot_ly(data = prostate, x = ~factor(svi), y = ~lpsa, 
                    type = "box", name="lpsa", width = 900) %>%
  layout(xaxis = list(title = "svi"), yaxis = list(title = "lpsa"))

fig <- subplot(boxplot1, boxplot2, boxplot3, nrows = 1, titleX = TRUE, titleY = TRUE) %>%
  layout(title="Distribución de los predictores según svi")
fig
```

Como podemos observar, la distribución de `svi` respecto al resto de variables es muy diferenciable. A simple vista podemos ver una clara relación entre valores negativos del `svi` con los tres predictores propuestos.

Vamos ahora a calcular el modelo. Dado que la variable `svi` es binaria (1 o 0), especificamos el parámetro `family=binomial`:
```{r}
# Nos aseguramos de que la variable svi esté en formato factor, por si se ha cambiado en algun apartado anterior
prostate$svi <- as.factor(prostate$svi)

modelo_logistico <- glm(svi ~ lcavol + lcp + lpsa, data = prostate, family = binomial)
# Creamos una lista con los coeficientes preparados para el display
coefs = sapply(modelo_logistico$coefficients, function(x) round(x,2))
```

El modelo resulta en la siguiente función:

$$
\text{svi} = `r coefs[1]` + `r coefs[2]` \cdot \text{lcavol} + `r coefs[3]` \cdot \text{lcp}+ `r coefs[4]` \cdot \text{lpsa}
$$

## 8.2 Conclusiones sobre el modelo

Veamos el `summary` del modelo:

```{r}
summary(modelo_logistico)
```

Lo primero que salta a la vista es que el error estimado para la variable `lcavol` es unas 6 veces mayor que el coeficiente en sí, lo cual nos indica que podría ser que la significancia de este predictor sea baja. Además, recordemos del **apartado 5** que `lcavol` guardaba una alta correlación con las otras variables, `lcp` y `lpsa`. Esto se nos confirma al observar el $p-$valor de esta variable, de aproximadamente $0.87$, lo cual nos indica que este predictor no es estadísticamente significativo. 

Vamos a realizar una predicción de nuestro set original usando el modelo:
```{r}
threshold = 0.5
log.probs = predict(modelo_logistico, type = "response")
log.preds = rep(0, length(log.probs))
log.preds[log.probs > threshold] = 1
```

Veamos la tabla de valores originales:
```{r}
table(log.preds)
```

Y a continuación la tabla de valores obtenidos:
```{r}
table(prostate$svi)
```

Como se puede observar, el modelo solo falló en 3 ocasiones, dando un $0$ cuando debería ser $1$, lo cual indica que el modelo hizo un buen trabajo al predecir el conjunto de entrenamiento, con un `r round((79+18-3)/(79+18)*100)`% de precisión. Es importante hacer notar que el dataset está cargado hacia `svi = 0`, es decir, hay muchas más observaciones en las que `svi = 1`, lo cual en general no es bueno para la generación de modelos, y en particular puede ser que haya influenciado a esta tendencia al 0 que muestra nuestro modelo. 

De cualquier forma, en esta ocasión no hemos dividido el dataset en *training* y *testing*, por lo que estamos examinando el rendimiento de nuestro modelo en el mismo conjunto de datos con el que lo hemos entrenado. En este caso habría que tener en cuenta dos cosas:

1. La precisión del modelo no es representativa ya que, como decimos, no estamos utilizando un conjunto de testeo.
2. Al analizar la precisión y obtener un resultado tan bueno, deberíamos tener cuidado con el **overfitting**.

Por estos motivos, se concluye que, a priori, no se puede decir que sea un buen modelo sin realizar primero pruebas de su rendimiento con un set de testeo, además de que mejoraría en gran medida si se entrenara con un set más equilibrado.

### Modelo alternativo

Vistos los resultados obtenidos con el modelo anterior, nos preguntamos si podríamos mejorarlo más. Y es que, como explicamos, el predictor `lcavol` tiene una baja significancia, esto junto con la alta correlación con el resto de predictores del modelo nos lleva a pensar que podemos mejorar el modelo al eliminar este predictor.

Vamos a proceder de forma similar al modelo original:
```{r}
modelo_logistico.nolcavol <- glm(svi ~ lcp + lpsa, data = prostate, family = binomial)
# Creamos una lista con los coeficientes preparados para el display
coefs = sapply(modelo_logistico.nolcavol$coefficients, function(x) round(x,2))
summary(modelo_logistico.nolcavol)
```

El modelo resulta en la siguiente función:
$$
\text{svi} = `r coefs[1]` + `r coefs[2]` \cdot \text{lcp}+ `r coefs[3]` \cdot \text{lpsa}
$$

Además, los $p-$values nos muestran que todas las variables utilizadas, `lcp` y `lpsa` en este caso, son significantes. Vamos a comparar los resultados obtenidos 

Vamos a realizar una predicción de nuestro set original usando el modelo:
```{r}
threshold = 0.5
log.probs = predict(modelo_logistico.nolcavol, type = "response")
log.preds = rep(0, length(log.probs))
log.preds[log.probs > threshold] = 1
```

Veamos la tabla de valores originales:
```{r}
table(log.preds)
```

Y a continuación la tabla de valores obtenidos:
```{r}
table(prostate$svi)
```

Observamos como hemos conseguido predecir un valor más utilizando este modelo, que se desprende de la variable `lcavol`, llegando así al `r round((79+18-2)/(79+18)*100)`% de precisión. Sin embargo, como también mencionamos con el modelo anterior, debemos ser críticos con este modelo y prestar especialmente atención al posible **overfitting** y al desbalance del conjunto de datos.

## 8.3 Predicción con el modelo

Vamos a utilizar nuestro modelo para predecir un valor de `svi` para valores concretos de los predictores `lcavol`, `lcp` y `lpsa`. Utilizamos los siguientes valores:

- `lcavol` = 2.8269,
- `lcp` = 1.843,
- `lpsa` = 3.285

Primero creamos el `data.frame` que contenga estos datos:
```{r}
pred_data = data.frame(
  lcavol = 2.8269,
  lcp = 1.843,
  lpsa = 3.285
)
```

Realizamos la predicción usando nuestro modelo:
```{r}
prob <- predict(modelo_logistico, newdata = pred_data, type = "response")
```

Obtenemos una probabilidad del `r round(prob*100,1)`% de que `svi`$=1$ con estos valores de los predictores. 

### Predicción usando el modelo alternativo

Vamos a realizar la misma predicción, pero esta vez utilizando el modelo propuesto en el subapartado anterior, que se desprendía del predictor `lcavol` por su baja significancia estadística. Dado que los datos ya están preparados, podemos calcular la probabilidad con un solo comando (se ignora el valor de `lcavol` al realizar la predicción):

```{r}
prob.nolcavol <- predict(modelo_logistico.nolcavol, newdata = pred_data, type = "response")
```

Obtenemos una probabilidad del `r round(prob.nolcavol*100,1)`% de que `svi`$=1$ con estos valores de los predictores. Podemos observar que la probabilidad es prácticamente la misma que con el modelo anterior.

# 9. PCA-PCR

## 9.1 Biplot y `summary` del modelo

Primero, creamos nuestro modelo de componentes principales **PCA**:
```{r}
# Creamos el subset para el PCA
pca.data <- prostate[,c("lcavol", "lweight", "lbph", "lcp")]
# Añadimos el svi a parte, ya que al ser categorico hay que reconvertirlo a numerico
pca.data["svi"] <- as.numeric(prostate$svi)
# Realizamos el PCA
pca.model <- prcomp(pca.data, scale. = TRUE, center = TRUE)
pca.var <- pca.model$sdev ^2
pca.pve <- pca.var / sum(pca.var)
```

### `summary` del modelo

```{r}
summary(pca.model)
```

Podemos observar en la primera fila del `summary` de nuestro modelo, *Standard deviation*, la cantidad de variabilidad que cada componente principal captura de los datos. Como se puede observar, los valores disminuyen con cada componente, siendo la mayor **PC1** y la menor **PC5**.

### Desviación estándar de la primera componente principal

La desviación estándar en un análisis de componentes principales representa la cantidad de variabilidad de los datos que es capturada por esta componente, lo cual es una medida clave en el PCA, ya que su objetivo es el de identificar las direcciones, o componentes principales, que mejor capturen la varianza del conjunto de datos.

El primer componente principal `PC1` es, por definición, aquel que caputra la mayor varianza posible, maximizando la varianza en las observaciones proyectadas en esa dirección, con lo cual una gran parte de las diferencias entre observaciones se ven reflejadas a lo largo de esta dirección.

Por tanto, se puede decir que la mayor desviación estándar de esta primera componente es consecuencia del objetivo de un PCA.

### Componentes principales

Las siguientes filas del `summary` ahondan en esta descripción, proporcionándonos valores porcentuales y acumulativos de la descripción de la varianza explicada por cada componente. Representemos estos valores gráficamente:
```{r}
par(mfrow = c(1, 2))
plot(pca.pve, xlab = "Principal Component",
  ylab = "Proportion of Variance Explained", ylim = c(0, 1),
  type = "b")
plot(cumsum(pca.pve), xlab = "Principal Component",
  ylab = "Cumulative Proportion of Variance Explained",
  ylim = c(0, 1), type = "b")
```

Como podemos observar, un $75%$ de la varianza de los datos originales se podrían explicar usando tan solo las dos primeras componentes principales, lo cual sugiere que estas dos componentes podrían ser suficientes para representar los datos, con una pequeña pérdida de información sobre la variación de los datos, llegando al $85%$ en caso de utilizar **tres componentes**. 

### `biplot` con las dos primeras componentes principales

Realicemos un **biplot** para ver cómo se expresarían nuestros datos en función de esas dos componentes.

```{r}
biplot(pca.model, scale = 0)
```

En este biplot podemos observar los puntos que representan a los pacientes en el espacio definido por las dos primeras componentes principales (`PC1` y `PC2`). 

Podemos observar diferentes agrupaciones de puntos, lo que indicaría que estas observaciones tendrían valores similares también en las variables originales.

Además, gracias a las flechas que indican la composición de las componentes principales según las variables originales, podemos observar que, al haber un ángulo pequeño entre ellas, las variables `lcavol`, `lcp` y `svi` están fuertemente correlacionadas, como se corrobora en apartados anteriores. Lo mismo ocurriría con `lbph` y `lweight` aunque en menor medida, al haber un ángulo de separación mayor entre ellas. 

## 9.2 PCR

Iniciamos con una carga de las librerias necesarias

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


### Análisis exploratorio de datos y correlación
Iniciamos echando una pequeña vista a la correlación existente entre los datos. Se ha eliminado la diagonal para mejorar la visualización.
```{r}
df <- prostate
df$svi <- as.numeric(df$svi)

corrplot(cor(df[, c("lcavol", "lweight", "lbph", "lcp", "svi")]), 
         type = "upper", 
         tl.cex = 0.7, 
         diag = FALSE,
         addCoef.col = "black") 
```

Si hubiéramos usado la función `findCorrelation`, no se vería ninguna hasta modificar el valor por defecto de `cutoff = 0.95`. 
```{r}
findCorrelation(cor(df[, c("lcavol", "lweight", "lbph", "lcp", "svi")]), 
                names = TRUE,
                cutoff = 0.44)
```
## Ajuste de la Regresión con Componentes Principales (PCR)
```{r}
suppressPackageStartupMessages({
  library(pls)
})
```

Para abordar la colinealidad, se ajusta un modelo de PCR y uso validación cruzada para evaluar el error de predicción. 

Se incluye `scale = TRUE` para normalizar los predictores y `validation = "CV"` para usar validación cruzada en la evaluación:

```{r}
set.seed(42)
pcr.fit = pcr(lpsa ~ lcavol + lweight + lbph + lcp + svi, 
              data = df, 
              scale = TRUE, 
              validation = "CV")
```


### Visualización del MSE y selección de componentes principales

Para decidir el número óptimo de componentes principales, se traza las puntuaciones de la validación cruzada utilizando `validationplot()`. Visualmente ya se ve que  el menor error de validación cruzada es con M = 5. 

Además se ve que el aumento de componentes favorece en los resultados. Haciendo uso completo de los pedidos en el enunciado. Aunque de forma relativa el uso de 1 a 4 componentes podría ser equivalente si quisiéramos usar menos.

```{r}
validationplot(pcr.fit, val.type = "MSEP")
```

Fijándonos en el `summary()` podemos ver mismamente que en `M = 5` se captura el 100% de la varianza en los predictores.

```{r}
summary(pcr.fit)
```

O a través de la función `which.min()` vemos el error más pequeño:
```{r}
which.min( MSEP( pcr.fit )$val[1,1, ] ) - 1
```
### Comparación con Regresión lineal múltiple

```{r}
# Predicciones
pcr.pred_full <- predict(pcr.fit, df, ncomp = 5)
multi_pred_full <- predict(modelo_lpsa, newdata = df)



# SSE para PCR
sse_pcr_full <- sum((pcr.pred_full - df$lpsa)^2)

# SSE para Regresión Lineal Múltiple
sse_multi_full <- sum((multi_pred_full - df$lpsa)^2)


n <- nrow(df)
p_pcr <- 5  
p_multi <- length(coef(modelo_lpsa)) 


# Calcular el RSE para PCR
rse_pcr_full <- sqrt(sse_pcr_full / (n - p_pcr))
# Calcular el RSE para Regresión Lineal Múltiple
rse_multi_full <- sqrt(sse_multi_full / (n - p_multi))




# R2 pcr
ss_total_full <- sum((df$lpsa - mean(df$lpsa))^2)
ss_residual_pcr_full <- sum((pcr.pred_full - df$lpsa)^2)
r_squared_pcr_full <- 1 - (ss_residual_pcr_full / ss_total_full)
adj_r_squared_pcr_full <- 1 - ((1 - r_squared_pcr_full) * (nrow(df) - 1) / (nrow(df) - 2 - 1))
# R2 multi
adj_r_squared_multi_full <- summary(modelo_lpsa)$adj.r.squared


cat("PCR: regresión de componentes principales\n",
    "MULTI: Regresión Lineal Múltiple\n\n",
    "############ RSE ############\n\n",
    "RSE (PCR):", rse_pcr_full, "\n",
    "RSE (MULTI):", rse_multi_full, "\n\n",
    "############ R² ############\n\n",
    "R² ajustado (PCR):", adj_r_squared_pcr_full, "\n",
    "R² ajustado (MULTI):", adj_r_squared_multi_full, "\n")
```
**RSE (Error Estándar de la Regresión)**:

Ambos modelos, **PCR** y **Regresión Lineal Múltiple (MULTI)**, tienen un RSE muy similar:

- **RSE (PCR)**: `r round(rse_pcr_full, 4)`  
- **RSE (MULTI)**: `r round(rse_multi_full, 4)`

Esto indica que ambos modelos tienen un error estándar de la regresión comparable, lo que sugiere que la dispersión de los errores de las predicciones es casi igual en ambos casos.

**R² ajustado**:

- El modelo de **PCR** tiene un **R² ajustado de 0.6368**, lo que significa que explica el **63.68%** de la variabilidad de los datos.
- El modelo de **Regresión Lineal Múltiple (MULTI)** tiene un **R² ajustado de 0.6248**, lo que implica que explica el **62.48%** de la variabilidad de los datos.

**Conclusión**:

Aunque ambos modelos presentan **RSE similares**, lo que refleja una dispersión de errores comparable, el modelo **PCR** tiene una ligera ventaja en cuanto a **R² ajustado**, indicando que **PCR** explica mejor la variabilidad de los datos. Este pequeño beneficio en el ajuste hace que **PCR** sea el modelo preferido en este caso, ya que demuestra un mayor rendimiento al ajustar los datos.

Sin embargo, la diferencia no es tan significativa, por lo que dependiendo del contexto y los requisitos del análisis, **Regresión Lineal Múltiple (MULTI)** también puede ser una opción válida si se busca simplicidad o interpretabilidad adicional.