---
title: "Partial Least Squares Regression (PLS)"
author: "Sonia Tarazona"
date: "`r [Link]()`"
output:
pdf_document:
toc: true
toc_depth: '2'
html_document:
toc: true
number_sections: true
toc_depth: 2
toc_float:
collapsed: false
smooth_scroll: true
---
<style>
.header-section-number::after { content: "."; }
</style>
```{r setup, warning=FALSE, message=FALSE, cache=FALSE}
library(knitr)
knitr::opts_chunk$set(echo = TRUE)
library(caret)
# [Link]("devtools")
# devtools::install_github("BiostatOmics/PLSandO")
# pak::pak("BiostatOmics/PLSandO")
library(PLSandO)
library(gridExtra)
```
# Preparación de datos para modelos PLS
**NOTA IMPORTANTE**
Al igual que en PCA, la librería *PLSandO* nos permite utilizar la función **Preparing()**
para procesar nuestros datos antes de aplicar cualquier modelo. Este preprocesamiento
incluye la conversión de factores a variables binarias 0-1 que puedan ser incluidas como
predictoras en los modelos PLS o PLS-DA, la eliminación de variables/observaciones con
exceso de valores faltantes o el filtrado de variables con baja o nula variabilidad.
Si no se aplica la función **Preparing()**, el propio modelo PCA, PLS o PLS-DA aplicará
automáticamente el preprocesado con los parámetros por defecto de dicha función.
Es importante tener en cuenta que los modelos PLS de *PLSandO* se ajustan mediante el
algoritmo NIPALS, que admite valores faltantes.
# PLS2: Aceites
## Lectura de datos
Cargamos los datos de ejemplo del aceite de oliva, que nos servirán de ejemplo para
realizar tanto un PLS2 como un PLS-DA. La matriz **X** contiene variables químicas y la
matriz respuesta **Y** contiene variables sensoriales. Además, el vector *proced* contiene
la procedencia de cada aceite (S=España, I=Italia, G=Grecia).
```{r datosAceite}
load("archivos de datos/[Link]", verbose = TRUE)
```
Para el modelo PLS2, utilizaremos la matrix **X** como matriz predictora, y la **Y** como
matriz respuesta.
## Estimación del modelo y optimización del número de componentes
Escalamos tanto la matriz **Y** como la **X**, ya que las variables están medidas en
distintas unidades. Lo haremos desde la propia función **pls()**, teniendo en cuenta que las
opciones seleccionadas en los argumentos correspondientes centran y escalan cada matriz.
No dividiremos los datos en entrenamiento y test porque tenemos muy pocas observaciones
($n =$ `r nrow(X)`).
Estimaremos el número de componentes óptimo mediante validación cruzada. En este
caso, al tener un número tan reducido de observaciones, optaremos por el procedimiento
"leave-one-out" (LOO).
```{r selComps, message = FALSE, warning=FALSE, [Link]=3.5, [Link]=3}
mypls = pls(X, Y, cvFolds = nrow(X), rep = 1, perm = 20,
scaling = "standard", scalingY = "standard",
algo = "nipals", train = 1, alpha = 0.05, parallel = FALSE)
kable(mypls$summary, digits = 3)
plsPlot(mypls, type = 'R2vsQ2')
```
El criterio de la función **pls()** para seleccionar el número óptimo de componentes se basa
en los valores de $R^2$ y de $Q^2$. Se extraen componentes mientras se cumplan las dos
condiciones siguientes: (1) la $R^2$ de la nueva componente supere en 0.01 la $R^2$ de la
anterior y (2) la diferencia entre la $Q^2$ de la nueva componente y la de la componente
anterior sea positiva. De acuerdo con este criterio, el número óptimo de componentes es `r
mypls$ncomp`.
Confirmamos en el gráfico de $Q^2$ y $R^2$ que, efectivamente, `r mypls$ncomp` es un
número de componentes adecuado ya que a partir de la tercera componente baja bastante
el valor de $Q^2$.
Podemos observar en el resumen mostrado que el modelo con `r mypls$ncomp`
componentes explica un `r round(mypls$summary[2,"cumR2Y"],2)*100`% de la variabilidad
de las variables en **Y** y un `r round(mypls$summary[2,"cumR2X"],2)*100`% de la
variabilidad de las variables en **X**, siendo la capacidad predictiva de $Q^2 =$ `r
round(mypls$summary[2,"cumQ2"],2)`. Además nos indica el valor del RMSE, calculado
como la suma de los RMSE de cada variable en $Y$. En este caso sirve para confirmar que
baja al añadir la segunda componente pero para evaluar mejor la bondad de las
predicciones del modelo, calcularemos más adelante el RMSE por variable y otras métricas
de error.
## Validación del modelo PLS2
### Detección de casos anómalos moderados con la SCR (distancia al modelo)
En los siguientes gráficos representamos la Suma de Cuadrados Residual (SCR) con
límites de confianza al 95% y al 99%.
```{r SCR, [Link]=7, [Link]=3.5}
par(mfrow = c(1,2))
anomalos = plsOutliers(mypls, method = "RSS", conf = 0.95)
anomalos99 = plsOutliers(mypls, method = "RSS", conf = 0.99)
```
**Conclusión:** En este caso, hay `r length(anomalos$ModerateOutliers)` aceite que se sale
fuera del límite del 95%. Esta observación mal explicada por el modelo corresponde al
aceite `r anomalos$ModerateOutliers`. No la excluiremos, puesto que ni siquiera excede el
límite del 99%. No obstante, generaremos el gráfico de contribuciones para esta
observación, de modo que nos ayude a entender el por qué de su anomalía:
```{r, [Link]=3, [Link]=3}
plsOutlierContrib(mypls, outliers = anomalos)
```
El aceite G4 tiene un valor anormalmente bajo en la variable *DK* y alto en *K270*.
----------
*EJERCICIO 1*
*Recupera los valores originales del aceite G4 en las variables DK y K270 y discute si son
errores de medida o si son valores que, aunque raros, es posible que se hayan dado en la
realidad. En caso de que sean errores de medida, ¿qué harías?*
----------
### Detección de anómalos severos con T2-Hotelling
Los siguientes gráficos para detectar anómalos severos muestran los valores de la $T^2$
de Hotelling para cada observación, con los límites de confianza al 95% y al 99%,
respectivamente.
```{r T2b, [Link]=7, [Link]=3.5}
par(mfrow = c(1,2))
severos = plsOutliers(mypls, method = "T2", conf = 0.95)
severos99 = plsOutliers(mypls, method = "T2", conf = 0.99)
```
En este caso, no hay ningún aceite que se salga fuera de estos límites, por lo que no hay
ningún anómalo severo.
### Relación lineal entre scores
Comprobaremos a continuación la relación de linealidad entre los scores de **X** y de **Y**
para cada componente mediante los siguientes gráficos, que incluyen también el coeficiente
de correlación lineal de Pearson.
```{r tu, [Link]=7, [Link]=3}
p1 = plsPlot(mypls, type = "linearity", comp = 1)
p2 = plsPlot(mypls, type = "linearity", comp = 2)
[Link](p1, p2, nrow = 1)
```
----------
*EJERCICIO 2*
*¿Podemos asumir que se cumple el supuesto de linealidad del modelo PLS?*
----------
### Cálculo de $R^2$ por componente y variable
Analizaremos con los siguientes gráficos el valor de $R^2$ variable a variable y por
componente, tanto para el espacio $X$ como para el espacio $Y$:
```{r R2, [Link]=7, [Link]=3.5}
p1 = plsPlot(mypls, type = "R2X")
p2 = plsPlot(mypls, type = "R2Y")
[Link](p1, p2, nrow = 1)
```
- En estos gráficos, estamos representando la proporción de variabilidad explicada por cada
componente de cada una de las variables tanto del espacio X como del espacio Y.
- En el espacio de las X, observamos que variables como *K270* o *K232* están
mayormente explicadas por la componente 1, al contrario que variables como *Acidity*.
- En el espacio de las Y, podemos observar que variables como *transp* o *glossy* están
mayormente explicadas por la componente 1. Sabemos que estarán muy correlacionadas.
### Validez del modelo
El siguiente gráfico de validación del modelo PLS se basa en la técnica de permutación. Se
permutan las filas de **Y** $perm=20$ veces (ver argumentos dados a la función **pls()**),
mientras **X** se deja invariable y se obtiene un modelo PLS para cada **Y** permutada
($Y_{perm}$). En cada modelo "permutado", se calculan los valores de la $R^2$ y de la
$Q^2$ y se representan con puntos morados y verdes, respectivamente. Las líneas
horizontales representan los valores de $R^2$ y $Q^2$ del modelo sin permutar **Y**. En el
eje X se representa la similitud (coeficiente de correlación de Pearson) que hay entre la
**Y** sin permutar y las $Y_{perm}$. En la parte superior del gráfico, los valores *pR2Y* y
*pQ2* muestran la proporción de permutaciones en las que los valores permutados mejoran
los obtenidos en el modelo original, sin permutar.
```{r, [Link]=4, [Link]=3.5}
plsPlot(mypls, type = "overfitting")
```
Los resultados muestran que tanto los valores de $R^2$ como de $Q^2$ del modelo PLS
obtenido con los datos sin permutar son mucho mejores que los obtenidos con los datos
permutados, por lo que podemos concluir que el modelo obtenido es mejor que un modelo
generado al azar y muestra evidencia de que no hay sobreajuste en el modelo PLS
ajustado.
## Interpretación del modelo PLS2
Una vez validado el modelo PLS2, pasamos a interpretarlo en profundidad. Para ello,
podemos analizar distintos gráficos, incluyendo los gráficos de scores en X y en Y, los
gráficos de loadings en X e Y (o de correlaciones), el gráfico de weights o el de coeficientes
de regresión.
```{r interpre, [Link]=7, [Link]=3}
# Score plots
p1 = plsPlot(mypls, type = "scoresX", comp = 1:2, colBy = proced)
p2 = plsPlot(mypls, type = "scoresY", comp = 1:2, colBy = proced)
[Link](p1, p2, nrow = 1)
# Loading plots
p1 = plsPlot(mypls, type = "loadingsX", comp = 1:2, colBy = "contrib")
p2 = plsPlot(mypls, type = "loadingsY", comp = 1:2, colBy = "contrib")
[Link](p1, p2, nrow = 1)
# Correlation plots
p1 = plsPlot(mypls, type = "corrX", comp = 1:2, colBy = "contrib")
p2 = plsPlot(mypls, type = "corrY", comp = 1:2, colBy = "contrib")
[Link](p1, p2, nrow = 1)
```
En primer lugar, vemos que, aunque no le proporcionamos al modelo la información sobre la
procedencia de cada aceite, las procedencias se separan bastante bien en las dos primeras
componentes, tanto en el espacio X como en el Y. La primera separa bien Italia y España,
mientras que la segunda ayuda a separar algunos aceites de Grecia. Así pues, a la vista del
gráfico de loadings o correlaciones en X, los aceites italianos tienden a tener más contenido
en *Peroxide* y *K232* que los españoles. Los griegos tienden a tener más acidez que el
resto.
----------
*EJERCICIO 3*
*¿Podrías caracterizar los aceites de cada país a partir del gráfico de loadings o
correlaciones en Y?*
----------
```{r interpre2, [Link]=3.5, [Link]=3.5}
# Weights
plsPlot(mypls, type = "weights", comp = 1:2, col = c("red3", "blue4"))
```
En segundo lugar, nos fijamos en el objetivo principal del PLS2, es decir, estudiar las
relaciones entre las variables en **X** y las variables en **Y**. Para ello, analizaremos
primero el gráfico de *weights*.
Entre otras cosas, podemos observar que valores más altos en *Peroxide* y *K232* se
corresponden con valores más altos en *brown* y *syrup*, y esto pasará especialmente en
aceites italianos. Los griegos (especialmente G1 y G4) tendrán valores más altos de
*Acidity* o *DK*, que a su vez están relacionados con valores más altos en *green*.
```{r interpre3, [Link]=6.5, [Link]=7.5}
# Coeficientes de regresión
plsPlot(mypls, type = "coef")
```
Por último, los gráficos de coeficientes de regresión nos permiten ver la relación de cada
variable en **X** con cada una de las variables en *Y*. Además, nos muestran la
significación estadística de dichos coeficientes obtenida mediante la técnica Jaccknife a
partir de la estrategia de validación cruzada utilizada en la estimación del modelo PLS, y las
variables están ordenadas de más a menos significativas. Un * significa que el p-valor
asociado al coeficiente de esa variable es menor que el nivel de significación del 5%.
----------
*EJERCICIO 4*
*A la vista de estos gráficos de coeficientes y de su significación estadística, ¿qué
conclusiones obtienes? ¿se confirman las observadas a partir del gráfico de weights?*
----------
## Selección de variables en PLS2
Aunque los gráficos de coeficientes anteriores nos permiten seleccionar las variables más
relevantes a partir de la significación estadística de los coeficientes obtenida mediante
Jackknife, *PLSandO* permite aplicar otras técnicas de selección de variables mediante la
función **plsVarSel()**: Jackknife mediante LOO, permutación, sMC (*Significance
Multivariate Correlation*) o VIP.
En el siguiente gráfico se muestra, por ejemplo, los resultados obtenidos mediante el
estadístico VIP, que muestra la importancia de cada variable en **X** para explicar **Y**:
```{r, [Link]=4, [Link]=3}
plsVarSel(mypls, type='VIP', threshold = 1)
```
Con los valores de VIP, confirmamos cuáles son las variables en **X** más importantes, en
general, para explicar **Y**: *Peroxide*, *K232* y *K270*, que son las 3 que sobrepasan el
valor mínimo de 1 "recomendado" para el VIP.
----------
*EJERCICIO 5*
*Utiliza la función plsVarSel() para estimar la significación estadística de los coeficientes de
regresión mediante la técnica de permutación con 1000 permutaciones. ¿Se obtienen los
mismos resultados que con la técnica Jackknife?*
----------
## Medidas del error de predicción
La función **pls()** nos proporciona algunas medidas del error como el RMSE pero de
forma global para todas las variables. Veamos cómo podemos calcular nuestras propias
medidas del error y hacerlo para cada variable en **Y** de forma independiente, ya que el
modelo podría estar explicando unas variables respuesta mejor que otras.
En primer lugar, con el siguiente código, podemos predecir los valores de la matriz
respuesta **Y** a partir del modelo PLS obtenido. Eso sí, como en este ejemplo no tenemos
datos *test*, solo podemos obtener las medidas del error sobre los datos de entrenamiento
que, como sabemos, no es lo más óptimo, ya que no serán medidas tan objetivas:
```{r pred1, [Link]=7, [Link]=3.5}
Ypred = predict(mypls, plot = FALSE)
residuos = Y-Ypred
myRMSE = sqrt(colMeans(residuos^2))
CVrmse = myRMSE/colMeans(Y)
par(mfrow = c(1,2))
barplot(myRMSE, las = 2, main = "RMSE")
barplot(CVrmse, las = 2, main = "CV-RMSE")
```
----------
*EJERCICIO 6*
*Comenta los resultados plasmados en los gráficos anteriores y discute la conveniencia de
usar RMSE o CV-RMSE para medir el error de predicción y comparar dicho error entre las
distintas variables respuesta.*
----------
Por último, los siguientes gráficos nos permiten ver de forma más detallada los valores
observados en cada variable de **Y** frente a los predichos por el modelo PLS, para saber
en qué casos se desvía más el modelo al hacer las predicciones.
```{r pred1b, [Link]=12, [Link]=8}
# Observados versus predichos
par (mfrow = c(2,3))
for (i in 1:ncol(Y)) {
plot(Y[,i], Ypred[,i], asp = 1, main = colnames(Y)[i],
xlab = "observado", ylab = "predicho", col = "white")
text(Y[,i], Ypred[,i], labels = rownames(Y))
abline(a=0, b=1, col = "red3", lwd = 2, lty = 2)
}
```
# PLS1: Cereales
## Lectura y preparación de datos
Cargamos los datos de ejemplo de cereales, ya procesados en prácticas anteriores. En este
caso, utilizaremos un modelo supervisado de regresión (PLS1) para explicar la variable
*rating* (**y**) en función de las variables relacionadas con el contenido nutricional de los
cereales (**X**).
```{r datosCereales}
load("[Link]", verbose = TRUE)
y = [Link]("rating" = cerePCA$rating)
X = subset(cerePCA, select = -rating)
```
## Estimación del modelo y del número de componentes
Escalamos la matriz **X**, ya que las variables nutricionales están medidas en distintas
unidades. Lo haremos desde la propia función **pls()**.
Dado que solo tenemos `r nrow(X)` observaciones, no dividiremos los datos en
entrenamiento y test, y estimaremos el número de componentes óptimo mediante validación
cruzada *k-fold*, con $k=10$, y 3 repeticiones. Además, paralelizaremos el cálculo del
modelo para que sea más eficiente.
```{r selCompsC, [Link]=3.5, [Link]=3, message=FALSE}
mypls = pls(x = X, y = y, ncomp = NULL, cvFolds = 10, rep = 3,
scaling = "standard", perm = 30, parallel = TRUE)
kable(mypls$summary, digits = 2)
```
Según los criterios de la librería *PLSandO*, el número óptimo de componentes sería `r
mypls$ncomp`.
------------
*EJERCICIO 7*
*Discute la bondad del modelo y medidas del error obtenidas.*
------------
## Validación del modelo PLS1
### Detección de casos atípicos
Aunque este análisis ya lo hicimos con el modelo PCA, lo repetiremos aquí para ilustrar
cómo se haría a partir del modelo PLS1.
En los siguientes gráficos representamos la Suma de Cuadrados Residual y la $T^2$ de
Hotelling con sus respectivos límites de confianza (99%).
```{r, [Link]=10, [Link]=5}
anomalos = plsOutliers(mypls, conf = 0.99)
anomalos$ModerateOutliers
anomalos$SevereOutliers
```
----------
*EJERCICIO 8*
*Discute los resultados obtenidos. ¿Consideras que hay algún caso atípico moderado o
severo que debiera excluirse? Genera los gráficos adecuados para entender a qué se
deben estas posibles anomalías.*
----------
### Relación lineal entre scores
```{r tuCer, [Link]=7, [Link]=3.5}
p1 = plsPlot(mypls, type = "linearity", comp = 1)
p2 = plsPlot(mypls, type = "linearity", comp = 2)
[Link](p1, p2, nrow = 1)
```
Tanto la relación lineal observada para en los gráficos para las dos componentes
(especialmente para la primera) como los respectivos coeficientes de correlación nos
permiten asumir la linealidad entre los scores de los espacios X e Y.
### Validez del modelo PLS1
```{r, [Link]=5, [Link]=3}
plsPlot(mypls, type = "overfitting")
```
----------
*EJERCICIO 9*
*¿Qué conclusiones extraerías del gráfico anterior?*
----------
## Interpretación del modelo PLS
Los siguientes gráficos nos servirán para interpretar el modelo PLS.
```{r, [Link]=8, [Link]=3.5}
# Score & loading plots
p1 = plsPlot(mypls, type = "scoresX", comp = 1:2,
colBy = [Link](shelf = estante), ellipses = FALSE)
p2 = plsPlot(mypls, type = "loadingsX", comp = 1:2, colBy = "contrib")
[Link](p1, p2, nrow = 1)
# Correlation in X and Weights
p1 = plsPlot(mypls, type = "corrX", comp = 1:2, colBy = "contrib")
p2 = plsPlot(mypls, type = "weights", comp = 1:2)
[Link](p1, p2, nrow = 1)
```
----------
*EJERCICIO 10*
*Interpreta los gráficos anteriores para extraer toda la información posible. ¿Qué variables
nutricionales son más determinantes en el valor del rating? ¿Por qué? ¿Son coherentes
estos resultados con los obtenidos anteriormente mediante modelos no supervisados?*
----------
## Selección de variables en el modelo PLS
El modelo PLS nos proporciona, por defecto, la significación estadística de los coeficientes
de regresión para cada variable predictora mediante el procedimiento *jackknife*:
```{r, [Link]=5, [Link]=3.5}
plsPlot(mypls, type = "coef", comp = 1:2)
```
Podemos también utilizar la función **plsVarSel()**, que nos permite estimar dicha
significación mediante otros procedimientos como, por ejemplo, la técnica de permutación, o
cambiar el nivel de significación:
```{r, [Link]=5, [Link]=3}
# Jackknife con LOO
plsVarSel(mypls, type = 'Jack', threshold = 0.05)
# Permutación
plsVarSel(mypls, type = 'Perm', rep = 1000, threshold = 0.05)
# VIP
plsVarSel(mypls, type = 'VIP', threshold = 1)
# Significance Multivariate Correlation
plsVarSel(mypls, type = 'sMC', threshold = 0.05)
```
----------
*EJERCICIO 11*
*A partir de los resultados anteriores, ¿cuáles crees que son las variables más relevantes
para predecir el rating de los cereales. Obtén un nuevo modelo PLS solo con las variables
que consideres que deban ser seleccionadas.*
----------
## Medidas del error de predicción en PLS1
----------
*EJERCICIO 12*
*Calcula las medidas del error que consideres más apropiadas para el modelo con todas las
variables y para el modelo con las variables seleccionadas obtenido en el ejercicio anterior.
Compara las medidas del error entre ambos modelos y discute los resultados.*
----------
# PLS-DA (2 clases): Cáncer
## Lectura y preparación de datos
Utilizaremos en esta sección los datos de cáncer de mama analizados en la Práctica 5 y los
dividiremos en datos de entrenamiento (80%) y test (20%), pero esta vez lo podemos hacer
con el propio paquete *PLSandO*.
```{r tumores, message=FALSE}
cancer = [Link]("[Link]", sep = ";", [Link] = 1)
nombres = c("rad", "tex", "per", "are", "smoo", "comp",
"conc", "cp", "sym", "frd")
colnames(cancer)[-1] = c(paste0(nombres,"_m"), paste0(nombres,"_se"),
paste0(nombres,"_p"))
cancer$diagnosis = factor(cancer$diagnosis)
y = [Link](diag = cancer$diagnosis)
X = subset(cancer, select = -diagnosis)
```
## Estimación del modelo PLS-DA y del número de componentes
En este caso, escalaremos la matriz **X** en la propia función *pls* y aplicaremos una
validación cruzada k-fold repetida para estimar el número de componentes.
```{r plsDAcan, message=FALSE, [Link]=12, [Link]=4}
[Link](100)
myplsda = PLSandO::plsda(x = X, y = y, ncomp = NULL, cvFolds = 10, rep = 3,
perm = 20, scaling = "standard", train = 0.8)
plsdaPlot(myplsda, type = 'ncomp')
```
El modelo estima en `r myplsda$ncomp` el número óptimo de componentes, aplicando en el
caso del PLS-DA la búsqueda del primer codo en la métrica *Balanced Error Rate* (BER).
Sin embargo, el F1-score y el criterio utilizado en PLS (regla de la $R^2$ y $Q^2$) sugieren
3 componentes. Dado que la variable respuesta en PLS-DA es categórica, este último
criterio no suele ser tan apropiado. Puede ser aconsejable representar gráficos de scores
para las 3 componente y asegurarnos que con solo dos podemos separar perfectamente los
tumores benignos y malignos.
La tabla siguiente muestra las medidas del error mostradas en los gráficos y alguna otra,
tanto para los datos de entrenamiento como también para los datos test. Todas las métricas
tienen valores bastante buenos.
```{r}
resumen = myplsda$summary
resumen[,-1] = round(resumen[,-1], 4)
kable(t(resumen))
```
## Validación del modelo PLS-DA
----------
*EJERCICIO 13*
*Estudia la existencia de valores anómalos mediante la $T^2$ de Hotelling y la SCR.
¿Consideras necesario excluir alguna observación?*
----------
*EJERCICIO 14*
*Valora la existencia de sobreajuste en el modelo obtenido.*
----------
## Interpretación del modelo PLS-DA
Aunque como hemos visto anteriormente se pueden obtener más gráficos, nos centraremos
en los gráficos de *scores* y *weights* para interpretar este modelo. Dado que el número de
variables predictoras es muy elevado, mostramos solo el 30% más relevante. Este 30% de
variables corresponde a las variables con mayor contribución en esas componentes, y la
contribución se calcula a partir de los loadings.
```{r interprCan, [Link]=8, [Link]=4}
p1 = plsdaPlot(myplsda, type = 'scoresX', comp = 1:2)
p2 = plsdaPlot(myplsda, type = 'weights', comp = 1:2,
selVars = 0.3, col = c("red3", "blue4"))
[Link](p1, p2, nrow = 1)
plsVarSel(myplsda, type = "Perm", rep = 1000, selVars = 0.8)
```
----------
*EJERCICIO 15*
*Interpreta los gráficos anteriores. Evalúa los parámetros oportunos (VIP, coeficientes de
regresión,...) para identificar los predictores más discriminantes entre los tipos de tumor.
¿Tendrán estos predictores discriminantes un valor más alto para los tumores benignos o
para los malignos?*
----------
## Medidas del error en PLS-DA
En este ejemplo, podemos calcular las medidas del error tanto para los datos de
entrenamiento como para los datos test.
Para los datos de entrenamiento:
```{r prediDAtrain, message=FALSE, warning=FALSE}
mypred = predict(myplsda, plot = FALSE)
mypred = factor(mypred$Class$diag)
myobs = myplsda$input$Y$diag
caret::confusionMatrix(mypred, myobs, positive = "M")
```
Para los datos test (o para nuevas observaciones futuras), podemos también obtener las
predicciones a partir del modelo PLS-DA ajustado. En este caso, es importante verificar que
las nuevas observaciones pertenecen a la misma población que la de los datos de
entrenamiento. Si no fuera así, no tendría sentido hacer predicciones sobre ellas. Para ello,
el gráfico que genera la función **predict()** proyecta las nuevas observaciones sobre el
espacio latente del modelo ajustado.
```{r prediDAtest, [Link]=4, [Link]=4}
mypred = predict(myplsda, new = myplsda$test$input$xTest, plot = TRUE)
```
Podemos observar en el gráfico que todas las observaciones de los datos test, excepto una
en la parte inferior derecha que se aleja en la segunda componente (menos relevante), se
proyectan sobre el espacio ocupado por los datos de entrenamiento y tiene sentido aplicar
sobre ellas el modelo ajustado para hacer predicciones.
Por tanto, veamos las medidas del error para los datos test. El objeto obtenido mediante la
función **predict()** contiene los siguientes elementos: `r names(mypred)`, que
corresponden a las predicciones numéricas para cada clase, las probabilidades de
pertenencia a cada clase (estimadas mediante el procedimiento *softmax*) y la clase
predicha por el modelo.
```{r prediDAtest2}
mypred = factor(mypred$Class$diag)
myobs = factor(myplsda$test$input$yTest$diag)
caret::confusionMatrix(mypred, myobs, positive = "M")
```
----------
*EJERCICIO 16*
*Discute las medidas del error obtenidas. ¿Mejoran dichas medidas de error si ajustamos un
modelo PLS-DA solo con los predictores más relevantes identificados en el Ejercicio 15?*
----------
*EJERCICIO 17*
*Obtén un gráfico con las curvas ROC para el modelo PLS-DA con todas las variables y
para el modelo de Análisis Discriminante obtenido en la práctica anterior y coméntalo.*
----------
# PLS-DA (3 clases): Aceites
Retomamos los datos de ejemplo del aceite de oliva que ya utilizamos al principio de esta
práctica. Nuestra variable respuesta será la procedencia del aceite (*proced*) y utilizaremos
como variables predictoras las dos matrices utilizadas en el modelo PLS (**X** e **Y**).
```{r datosAceiteDA, comment=FALSE}
load("archivos de datos/[Link]")
Xda = cbind(X, Y)
y = [Link]("proced" = proced)
```
## Estimación del modelo PLS-DA
Como disponemos de pocas observaciones, al igual que hicimos en el modelo PLS, no
dividiremos los datos en entrenamiento y test y utilizaremos una validación cruzada LOO
para estimar el número óptimo de componentes.
```{r plsDA, message=FALSE, warning=FALSE, [Link]=9, [Link]=3}
myplsda = plsda(x = Xda, y = y, ncomp = NULL, cvFolds = nrow(Xda), rep = 1,
perm = 40, algo = "nipals", train = 1, alpha = 0.05,
scaling = "standard", parallel = FALSE)
plsdaPlot(myplsda, type = 'ncomp')
```
A la vista de estos resultados, todas las métricas coinciden en que `r myplsda$ncomp`
componentes es una elección razonable para separar bien las procedencias del aceite.
## Validación del modelo
### Detección de casos atípicos
En los siguientes gráficos representamos la Suma de Cuadrados Residual y la $T^2$ de
Hotelling con sus respectivos límites de confianza (99%).
```{r, [Link]=7, [Link]=3.5}
anomalos = plsOutliers(myplsda, conf = 0.99)
```
No hay ningún aceite que se salga fuera de los límites del 99%, por lo que no hay ningún
anómalo moderado ni severo.
### Validez del modelo PLS-DA
```{r, [Link]=6, [Link]=4}
plsPlot(myplsda, type = "overfitting")
```
Como podemos observar, tanto los valores de $R^2$ como de $Q^2$ del modelo PLS-DA
obtenido con los datos sin permutar son mejores que los obtenidos con los datos
permutados, por lo que podemos concluir que el modelo obtenido es mejor que un modelo
generado al azar y no parece haber sobreajuste en el modelo PLS-DA ajustado.
## Interpretación del modelo
```{r interpr, message=FALSE}
p1 = plsdaPlot(myplsda, type = 'scoresX', comp = 1:2)
p2 = plsdaPlot(myplsda, type = 'weights', comp = 1:2, col = c("red3", "blue4"))
[Link](p1, p2, nrow = 1)
```
----------
*EJERCICIO 18*
*Interpreta los gráficos anteriores indicando las características de los aceites según su
procedencia.*
----------
## Selección de variables en el modelo PLS-DA
Al igual que en el modelo PLS, el modelo PLS-DA nos proporciona los valores de los
coeficientes de regresión y, por defecto, un intervalo de confianza para testar su
significación estadística mediante el procedimiento *jackknife*. Como tenemos 3 clases, se
obtienen los coeficientes de regresión para cada una de ellas:
```{r, [Link]=12, [Link]=8}
plsPlot(myplsda, type = "coef")
```
Análogamente, podemos utilizar la función **plsVarSel()** para estimar dicha significación
mediante otros procedimientos. En este caso, exploraremos, por ejemplo, el VIP:
```{r, [Link]=12, [Link]=4}
plsVarSel(myplsda, type = 'VIP', threshold = 1)
```
----------
*EJERCICIO 19*
*A la vista de todos los resultados anteriores, ¿qué variables consideras más relevantes
para predecir la procedencia del aceite?*
----------
## Medidas del error en PLS-DA
Obtendremos las medidas del error para los aceites con los que se ha entrenado el modelo,
ya que no disponemos de más datos.
```{r prediDA, message=FALSE, warning=FALSE}
mypred = predict(myplsda, plot = FALSE)
mypred = factor(mypred$Class$proced)
caret::confusionMatrix(mypred, proced)
```
----------
*EJERCICIO 20*
*Compara estos resultados con los obtenidos tras ajustar un modelo PLS-DA solo con las
variables seleccionadas en el Ejercicio 19.*
----------