Selección de predictores, regularización ridge, lasso, elastic net y reducción de dimensionalidad
Joaquín Amat Rodrigo
Diciembre 2016 (última actualización Noviembre 2020)
Más sobre ciencia de datos: cienciadedatos.net
- Regresión lineal simple
- Regresión lineal múltiple
- Validación de modelos predictivos: Cross-validation, OneLeaveOut, Bootstraping
- PCA
- Regularización Ridge, Lasso y Elastic Net con Python
Versión PDF: Github
Introducción
La regresión lineal múltiple es un método estadístico que trata de modelar la relación entre una variable continua y dos o más variables independientes mediante el ajuste de una ecuación lineal. Tres de las limitaciones que que aparecen en la práctica al tratar de emplear este tipo de modelos (ajustados por mínimos cuadrados ordinarios) son:
Se ven perjudicados por la incorporación de predictores correlacionados.
No realizan selección de predictores, todos los predictores se incorporan en el modelo aunque no aporten información relevante. Esto suele complicar la interpretación del modelo y reducir su capacidad predictiva. Existen otros modelos como random forest o gradient boosting que sí son capaces de seleccionar predictores.
No pueden ajustarse cuando el número de predictores es superior al número de observaciones.
Algunas de las estrategias que se pueden aplicar para atenuar el impacto de estos problemas son:
Subset selection: utilizar un proceso iterativo que vaya descartando los predictores menos relevantes.
Regularización Ridge, Lasso o Elastic Net: estos métodos fuerzan a que los coeficientes del modelo tiendan a cero, minimizando así el riesgo de overfitting, reduciendo varianza, atenuado el efecto de la correlación entre predictores y reduciendo la influencia en el modelo de los predictores menos relevantes.
Reducción de dimensionalidad: crean un número reducido de nuevos predictores (componentes) a partir de combinaciones lineales o no lineales de las variables originales y con ellas se ajusta el modelo.
A lo largo de este documento, se describen de forma progresiva los fundamentos teóricos y aspectos prácticos de cómo combinar regresión lineal múltiple con regularización, seleccion de predictores y reducción de dimensionalidad, y ejemplos de cómo crear este tipo de modelos en R.
Subset Selection
Los métodos conocidos como subset selection tienen la finalidad de identificar y seleccionar, de entre todos los predictores disponibles, aquellos que están más relacionados con la variable respuesta y así crear el mejor modelo. El esquema general de los métodos de subset selection consiste en:
Crear un conjunto de modelos candidatos (todos los posibles o un conjunto considerable de ellos), mediante diferentes combinaciones de los predictores disponibles.
Para cada posible tamaño de modelo (1 predictor, 2 predictores…) se selecciona el mejor basándose en el error de entrenamientor.
Los modelos finalistas de cada tamaño se comparan entre ellos para una métrica de validación (validación cruzada, \(C_p\), \(AIC\), \(BIC\) o \(R^2_{ajustado}\)) \(^{\text{Anexo 1}}\) y se identificar el mejor.
Dentro de esta estrategia se diferncian: best subset selection y stepwise selection (forward, backward e hybrid). Es importante tener en cuenta que, para un mismo conjunto de datos, no todos tienen por qué converger en un mismo modelo final.
Best subset selection
El proceso de best subset selection consiste en evaluar todos los posibles modelos que se pueden crear por combinación de los predictores disponibles. El algoritmo a seguir para \(k\) predictores es:
Se genera lo que se conoce como modelo nulo (\(M_0\)), que es el modelo sin ningún predictor.
Se generan todos los posibles modelos que contienen un único predictor y se selecciona el que tiene menor error de entrenamiento. Al modelo seleccionado se denomina (\(M_1\)).
Se repite el paso anterior para modelos con dos predictores y así sucesivamente hasta llegar al modelo con todos los predictores (\(M_k\)).
De entre los mejores modelos seleccionados para cada número de predictores (\(M_0\), \(M_1\), \(M_2\),…,\(M_k\)) se identifica el mejor modelo, esta vez empleando una métrica de validación (validación cruzada, \(C_p\), \(AIC\), \(BIC\) o \(R^2_{ajustado}\)). \(^{\text{Anexo 1}}\)
A pesar de que este método explora todas las posibilidades, tiene dos limitaciones fundamentales:
Rrequerimientos computacionales: Se requiere calcular \(2^p\) modelos distintos, lo que lo hace inviable para más de 40 predictores.
Problemas de overfitting. Al generarse tantos modelos, por simple azar se pueden encontrar buenos resultados. Por esta razón best subset selection no se ecominda si hay más de 10 predictores.
Forward stepwise selection
Forward stepwise selection es una alternativa computacionalmente más eficiente que best subset selection, en la que no se evalúan todas las posibles combinaciones de predictores sino solo un subconjunto. El proceso se inicia generando el modelo nulo (\(M_0\)) sin predictores. A continuación, se generan todos los posibles modelos que se pueden crear añadiendo un predictor al modelo nulo. De entre todos estos modelos con 1 predictor se selecciona el mejor basándose en el error de entrenamiento, al modelo elegido se denomina \(M_1\). Se repite el paso anterior, pero esta vez partiendo del último modelo seleccionado y así sucesivamente hasta llegar al modelo con todos los predictores. De entre los mejores modelos seleccionados para cada número de predictores (\(M_0\), \(M_1\), \(M_2\),…,\(M_k\)), se identifica el mejor, esta vez empleando una métrica de validación (validación cruzada, \(C_p\), \(AIC\), \(BIC\) o \(R^2_{ajustado}\)). \(^{\text{Anexo 1}}\)
Al crear modelos anidados, en los que el modelo \(k\) se construye a partir del modelo \(k-1\), el método forward stepwise selection no garantiza que se seleccione el mejor modelo de entre todos los posibles, ya que no se evalúan todas las posibles combinaciones. Sin embargo, suele llegar a modelos óptimos consiguiendo un buen rendimiento computacional y evitando el overfitting. Otra ventaja añadida es que, forward stepwise selection puede emplearse incluso cuando el número de predictores es mayor que el de observaciones.
Backward stepwise selection
El concepto es equivalente al de forward stepwise selection pero, en este caso, iniciando el proceso a partir del modelo que contiene todos los posibles predictores (full model \(M_k\)). En cada iteración, se entrenan todos los modelos que se pueden crear eliminando un único predictor y se selecciona el que tiene menor error de entrenamiento tiene. El proceso se repite hasta llegar al modelo nulo sin predictores (\(M_0\)). De entre los mejores modelos seleccionados para cada número de predictores (\(M_0\), \(M_1\), \(M_2\),…,\(M_k\)) se identifica el mejor, esta vez empleando una métrica de validación (validación cruzada, \(C_p\), \(AIC\), \(BIC\) o \(R^2_{ajustado}\)). \(^{\text{Anexo 1}}\) Backward stepwise selection permite evaluar cada variable en presencia de las otras, lo que es una ventaja frente a forward stepwise selection. Sin embargo, dado que el método se inicia con el modelo que contiene todos los predictores, si la regresión es por mínimos cuadrados, no se puede aplicar cuando el número de predictores es mayor que el número de observaciones.
Hibrid (double) Stepwise Selection
Este método se inicia al igual que el forward pero, tras cada nueva incorporación, se realiza un test de extracción de predictores no útiles (como en el backward). Este método se aproxima más al best subset selection pero sin caer en tantas limitaciones computacionales.
Regularización (Shrinkage)
Los métodos de subset selection descritos anteriormente emplean mínimos cuadrados ordinarios (OLS) para ajustar un modelo lineal que contiene únicamente un subconjunto de predictores. Otra alternativa, conocida como regularización o shrinkage, consiste en ajustar el modelo incluyendo todos los predictores pero aplicando una penalización que fuerce a que las estimaciones de los coeficientes de regresión tiendan a cero. Con esto se intenta evitar overfitting, reducir varianza, atenuar el efecto de la correlación entre predictores y minimizar la influencia en el modelo de los predictores menos relevantes. Por lo general, aplicando regularización se consigue modelos con mayor poder predictivo (generalización). Tres de los métodos de regularización más empleados son Ridge, Lasso y Elastic net.
Dado que estos métodos de regularización actúan sobre la magnitud de los coeficientes del modelo, todos deben de estár en la misma escala, por esta razón es necesario estandarizar o normalizar los predictores antes de entrenar el modelo. Ls métodos están especialmente indicados para situaciones en las que hay un mayor número de predictores que de observaciones.
Ridge
La regularización Ridge penaliza la suma de los coeficientes elevados al cuadrado \((||\beta||^2_2 = \sum_{j=1} ^p \beta^2_j)\). A esta penalización se le conoce como l2 y tiene el efecto de reducir de forma proporcional el valor de todos los coeficientes del modelo pero sin que estos lleguen a cero. El grado de penalización está controlado por el hiperparámetro \(\lambda\). Cuando \(\lambda = 0\), la penalización es nula y el resultado es equivalente al de un modelo lineal por mínimos cuadrados ordinarios (OLS). A medida que \(\lambda\) aumenta, mayor es la penalización y menor el valor de los predictores.
\[\sum^n_{i=1}(y_i - \beta_0 - \sum^p_{j=1} \beta_j x_{ij})^2 + \lambda \sum^p_{j=1} \beta_j^2 = \text{suma residuos cuadrados} + \lambda \sum^p_{j=1} \beta_j^2\]
La principal ventaja de aplicar ridge frente al ajuste por mínimos cuadrados ordinarios (OLS) es la reducción de varianza. Por lo general, en situaciones en las que la relación entre la variable respuesta y los predictores es aproximadamente lineal, las estimaciones por mínimos cuadrados tienen poco bias pero aún pueden sufrir alta varianza (pequeños cambios en los datos de entrenamiento tienen mucho impacto en el modelo resultante). Este problema se acentúa conforme el número de predictores introducido en el modelo se aproxima al número de observaciones de entrenamiento, llegando al punto en que, si \(p>n\), no es posible ajustar el modelo por mínimos cuadrados ordinarios. Empleando un valor adecuado de \(\lambda\), el método de ridge es capaz de reducir varianza sin apenas aumentar el bias, consiguiendo así un menor error total.
La desventaja del método ridge es que, el modelo final, incluye todos los predictores. Esto es así porque, si bien la penalización fuerza a que los coeficientes tiendan a cero, nunca llegan a ser exactamente cero (solo si \(\lambda=\infty\)). Este método consigue minimizar la influencia sobre el modelo de los predictores menos relacionados con la variable respuesta pero, en el modelo final, van a seguir apareciendo. Aunque esto no supone un problema para la precisión del modelo, sí lo es para su interpretación.
Lasso
La regularización Lasso penaliza la suma del valor absolutos de los coeficientes de regresión \((||\beta||_1 = \sum_{j=1} ^p |\beta_j|)\). A esta penalización se le conoce como l1 y tiene el efecto de forzar a que los coeficientes de los predictores tiendan a cero. Dado que un predictor con coeficiente de regresión cero no influye en el modelo, lasso consigue excluir los predictores menos relevantes. Al igual que en ridge, el grado de penalización está controlado por el hiperparámetro \(\lambda\). Cuando \(\lambda = 0\), el resultado es equivalente al de un modelo lineal por mínimos cuadrados ordinarios. A medida que \(\lambda\) aumenta, mayor es la penalización y más predictores quedan excluidos.
\[\sum^n_{i=1}(y_i - \beta_0 - \sum^p_{j=1} \beta_j x_{ij})^2 + \lambda \sum^p_{j=1} |\beta_j| = \text{suma residuos cuadrados} + \lambda \sum^p_{j=1} |\beta_j|\]
Comparación Ridge y Lasso
La principal diferencia práctica entre lasso y ridge es que el primero consigue que algunos coeficientes sean exactamente cero, por lo que realiza selección de predictores, mientras que el segundo no llega a excluir ninguno. Esto supone una ventaja notable de lasso en escenarios donde no todos los predictores son importantes para el modelo y se desea que los menos influyentes queden excluidos.
Por otro lado, cuando existen predictores altamente correlacionados (linealmente), ridge reduce la influencia de todos ellos a la vez y de forma proporcional, mientras que lasso tiende a seleccionar uno de ellos, dándole todo el peso y excluyendo al resto. En presencia de correlaciones, esta selección varía mucho con pequeñas perturbaciones (cambios en los datos de entrenamiento), por lo que, las soluciones de lasso, son muy inestables si los predictores están altamente correlacionados.
Para conseguir un equilibrio óptimo entre estas dos propiedades, se puede emplear lo que se conoce como penalización elastic net, que combina ambas estrategias.
Elastic net
Elastic net incluye una regularización que combina la penalización l1 y l2 \((\alpha \lambda ||\beta||_1 + \frac{1}{2}(1- \alpha)||\beta||^2_2)\). El grado en que influye cada una de las penalizaciones está controlado por el hiperparámetro \(\alpha\). Su valor está comprendido en el intervalo [0,1]. Cuando \(\alpha = 0\), se aplica ridge y cuando \(\alpha = 1\) se aplica lasso. La combinación de ambas penalizaciones suele dar lugar a buenos resultados. Una estrategia frecuentemente utilizada es asignarle casi todo el peso a la penalización l1 (\(\alpha\) muy próximo a 1) para conseguir seleccionar predictores y un poco a la l2 para dar cierta estabilidad en el caso de que algunos predictores estén correlacionados.
\[\frac{\sum^n_{i=1}(y_i - \beta_0 - \sum^p_{j=1} \beta_j x_{ij})^2}{2n} + \lambda (\alpha \sum^p_{j=1} |\beta_j| + \frac{1-\alpha}{2} + \sum^p_{j=1} \beta_j^2\]
Lambda Search
Encontrar el mejor modelo implica identificar el valor óptimo del hiperparámetro de regularización lambda (\(\lambda\)). Al tratarse de un hiperparámetro, no hay forma de saber de antemano qué valor es el adecuado. Una forma de lograrlo es emplear validación cruzada o generalized cross validation, esta última es una adaptación eficiente de leave-one-out cross validation disponible para la regulación Ridge.
Las principales implementaciones en R de estos métodos (glmnet y H2O) incorporan incorporan validación cruzada, sin embargo, hay situaciones para las que no son adecuadas. Por ejemplo, si los datos tienen un orden temporal (series temporales), la validación cruzada debe hacer un reparto ordenado de las observaciones, lo que se consigue fácilmente con otros paquetes como tidymodels o mlr3.
Aunque el valor óptimo de \(\lambda\) es aquel con el que se minimiza el error de validación cruzada, una práctica extendida es utilizar, en lugar de este, el mayor valor de \(\lambda\) que se aleja menos de una desviación típica del óptimo. De este modo, se consigue un modelo más sencillo (excluye más predictores) pero cuya capacidad predictiva es similar a la conseguida con el modelo más complejo.
Reducción de dimensionalidad
Los métodos anteriormente descritos controlan la varianza y overfitting, bien empleando únicamente un subconjunto de predictores, o bien haciendo que los coeficientes de regresión tiendan a cero. En ambos casos, se emplean las variables originales sin ser modificadas o como máximo habiendo sido estandarizadas.
Las técnicas conocidas como reducción de dimensionalidad crean un número reducido de nuevas variables (componentes) a partir de combinaciones lineales o no lineales de las variables originales, y con ellas se ajusta el modelo. De este modo, se consigue generar modelos con menor número de predictores pero que abarcan casi la misma información que la que aportan todas las variables originales. Existen diferentes aproximaciones para lograrlo, dos de las más utilizadas son: PCA y Partial Least Square (PLS).
Por lo general, ambos métodos tienen buenos resultados en aquellos casos en los que los predictores están altamente correlacionados. Cuando esto ocurre, con unas pocas componentes principales se puede capturar la mayor parte de la relación entre los predictores la variable respuesta. Es importante tener en cuenta que, aunque permiten generar modelos que contienen un número menor de predictores, no se trata de un método de selección de variables, ya que las nuevas variables (componentes) son combinaciones lineales de todos los predictores originales.
A continuación se describen de forma muy superficial los métodos de PCR y PLS. Para información más detallada sobre reducción de dimensionalidad consultar Análisis de Componentes Principales (Principal Component Analysis, PCA)
Principal Components Regression PCR
Principal Components Regression consiste en ajustar un modelo de regresión lineal mediante mínimos cuadrados empleando como predictores las componentes generadas a partir de un Principal Component Analysis (PCA). De esta forma, con un número reducido de componentes se puede explicar la mayor parte de la varianza. Algunas consideraciones a tener en cuenta son:
Cuando el número de componentes es igual al número de predictores originales, el resultado de Principal Components Regression es equivalente al de regresión por mínimos cuadrados ordinarios (OLS).
El número óptimo de componentes principales se puede elegir mediante validación cruzada.
Cuando se realiza PCR hay que estandarizar los predictores antes de calcular las PCA, de lo contrario, las variables que se miden en una escala mayor o las que presenten mayor varianza tendrán más peso. Si todos los predictores se miden con la misma escala y presentan la misma varianza, entonces no es necesaria la estandarización.
Partial least squares (PLS)
El método Partial Least Squares (PLS) es muy similar al PCR, en cuanto que ambos emplean como predictores las componentes principales de un PCA. La diferencia es que, mientras PCR ignora la variable respuesta Y para determinar las combinaciones lineales, PLS busca aquellas que, además de explicar la varianza observada, predicen Y lo mejor posible. Puede considerarse como una versión supervisada de PCR.
Ejemplo
El departamento de calidad de una empresa de alimentación se encarga de medir el contenido en grasa de la carne que comercializa. Este estudio se realiza mediante técnicas de analítica química, un proceso relativamente costoso en tiempo y recursos. Una alternativa que permitiría reducir costes y optimizar tiempo es emplear un espectrofotómetro (instrumento capaz de detectar la absorbancia que tiene un material a diferentes tipos de luz en función de sus características) e inferir el contenido en grasa a partir de sus medidas.
Antes de dar por válida esta nueva técnica, la empresa necesita comprobar qué margen de error tiene respecto al análisis químico. Para ello, se mide el espectro de absorbancia a 100 longitudes de onda en 215 muestras de carne, cuyo contenido en grasa se obtiene también por análisis químico, y se entrena un modelo con el objetivo de predecir el contenido en grasa a partir de los valores dados por el espectrofotómetro. Los datos meatspec empleados en este ejemplo se han obtenido del paquete faraway.
Paquetes
Los paquetes empleados en este ejemplo son:
# Datos
# ==============================================================================
library(faraway)
# Gráficos y tratamiento de datos
# ==============================================================================
library(tidyverse)
library(skimr)
#devtools::install_github("boxuancui/DataExplorer")
library(DataExplorer)
library(scales)
library(corrr)
# Modelado
# ==============================================================================
library(glmnet)
library(pls)Datos
El set de datos contiene 101 columnas. Las 100 primeras, nombradas como \(V1\) , …, \(V100\) recogen el valor de absorbancia para cada una de las 100 longitudes de onda analizadas (predictores), y la columna fat el contenido en grasa medido por técnicas químicas (variable respuesta).
# Datos
# ==============================================================================
data("meatspec")
datos <- meatspec
head(datos,3)Relación entre variables
El primer paso a la hora de establecer un modelo lineal múltiple es estudiar la relación que existe entre variables. Esta información es crítica a la hora de identificar cuáles pueden ser los mejores predictores para el modelo, y para detectar colinealidad entre predictores. A modo complementario, es recomendable representar la distribución de cada variable mediante histogramas.
# Correlación entre columnas numéricas
# ==============================================================================
df_correlaciones <- datos %>%
correlate(method = "pearson") %>%
stretch(remove.dups = TRUE)
df_correlaciones %>% mutate(r_abs = abs(r)) %>% arrange(desc(r_abs)) %>% head(5)Muchas de las variables están altamente correlacionadas (correlación absoluta > 0.8), lo que supone un problema a la hora de emplear modelos de regresión lineal.
Modelos
Se ajustan varios modelos lineales con y sin regularización, con el objetivo de identificar cuál de ellos es capaz de predecir mejor el contenido en grasa de la carne en función de las señales registradas por el espectrofotómetro.
Para poder evaluar la capacidad predictiva de cada modelo, se dividen las observaciones disponibles en dos grupos: uno de entrenamiento (70%) y otro de test (30%).
# División de los datos en train y test
# ==============================================================================
set.seed(1235)
id_train <- sample(1:nrow(datos), size = 0.7*nrow(datos), replace = FALSE)
datos_train <- datos[id_train, ]
datos_test <- datos[-id_train, ]Mínimos cuadrados (OLS)
En primer lugar, se ajusta un modelo de regresión lineal (OLS) incluyendo todas las longitudes de onda como predictores.
# Creación y entrenamiento del modelo
# ==============================================================================
modelo <- lm(fat ~ ., data = datos_train)
summary(modelo)##
## Call:
## lm(formula = fat ~ ., data = datos_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.82918 -0.38985 0.01879 0.43654 2.04621
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.885 3.185 2.475 0.01682 *
## V1 8992.864 5135.663 1.751 0.08619 .
## V2 -10649.370 9342.756 -1.140 0.25989
## V3 2876.109 14317.422 0.201 0.84162
## V4 2367.215 26770.078 0.088 0.92990
## V5 1335.792 35389.986 0.038 0.97004
## V6 -16468.795 29401.964 -0.560 0.57795
## V7 3254.574 17824.233 0.183 0.85587
## V8 -6107.216 11636.352 -0.525 0.60206
## V9 9222.489 8871.189 1.040 0.30363
## V10 10797.541 9583.731 1.127 0.26538
## V11 -17498.012 15182.337 -1.153 0.25470
## V12 36205.862 27881.944 1.299 0.20018
## V13 -35394.626 37620.319 -0.941 0.35140
## V14 23245.362 35357.095 0.657 0.51397
## V15 8464.127 26419.024 0.320 0.75004
## V16 -37353.236 16946.593 -2.204 0.03224 *
## V17 13688.144 10577.169 1.294 0.20169
## V18 5723.458 10395.573 0.551 0.58443
## V19 -6926.451 15432.974 -0.449 0.65555
## V20 32790.720 24364.409 1.346 0.18455
## V21 -82403.104 33862.845 -2.433 0.01865 *
## V22 88160.072 36171.668 2.437 0.01848 *
## V23 -40058.494 26845.339 -1.492 0.14206
## V24 4262.500 18180.180 0.234 0.81561
## V25 5858.536 11999.288 0.488 0.62756
## V26 -188.519 9113.771 -0.021 0.98358
## V27 -16091.298 10344.803 -1.555 0.12626
## V28 25189.662 14436.958 1.745 0.08729 .
## V29 -13930.381 17318.145 -0.804 0.42506
## V30 -1236.733 26545.480 -0.047 0.96303
## V31 -12272.314 30084.357 -0.408 0.68510
## V32 29011.792 26980.797 1.075 0.28752
## V33 -11818.513 21170.474 -0.558 0.57921
## V34 -5347.222 16151.217 -0.331 0.74200
## V35 -10547.717 13533.518 -0.779 0.43951
## V36 18920.496 11164.406 1.695 0.09648 .
## V37 -8377.894 10157.312 -0.825 0.41347
## V38 -8266.572 10876.386 -0.760 0.45087
## V39 28390.694 17148.511 1.656 0.10420
## V40 -27513.870 26702.791 -1.030 0.30789
## V41 19956.576 30947.110 0.645 0.52203
## V42 -6916.956 32724.987 -0.211 0.83348
## V43 -14498.503 29103.314 -0.498 0.62059
## V44 -3703.608 20350.036 -0.182 0.85634
## V45 34789.183 19266.416 1.806 0.07711 .
## V46 -16613.480 12909.347 -1.287 0.20416
## V47 -7777.822 5364.359 -1.450 0.15346
## V48 2196.605 6629.730 0.331 0.74181
## V49 4837.740 10209.503 0.474 0.63771
## V50 -3480.071 14069.281 -0.247 0.80567
## V51 -78.433 22816.541 -0.003 0.99727
## V52 8022.376 26344.945 0.305 0.76203
## V53 -29438.114 23122.364 -1.273 0.20897
## V54 51184.944 17339.904 2.952 0.00484 **
## V55 -47259.801 13783.140 -3.429 0.00124 **
## V56 25490.797 9004.096 2.831 0.00671 **
## V57 -8702.084 6652.006 -1.308 0.19691
## V58 -1545.866 6229.912 -0.248 0.80507
## V59 -234.443 7270.375 -0.032 0.97441
## V60 3159.078 6288.764 0.502 0.61768
## V61 2756.859 5530.132 0.499 0.62035
## V62 -1253.152 6072.163 -0.206 0.83735
## V63 11834.865 8209.067 1.442 0.15575
## V64 -20446.787 11554.133 -1.770 0.08301 .
## V65 4869.320 14888.725 0.327 0.74503
## V66 14359.138 17656.737 0.813 0.42002
## V67 -8161.971 17957.047 -0.455 0.65146
## V68 -20571.527 17588.322 -1.170 0.24781
## V69 21410.851 14938.916 1.433 0.15814
## V70 2546.365 11527.676 0.221 0.82609
## V71 -18824.174 12012.386 -1.567 0.12354
## V72 5484.490 12190.582 0.450 0.65477
## V73 16821.862 9147.796 1.839 0.07199 .
## V74 -10317.345 8521.342 -1.211 0.23179
## V75 -657.657 8727.037 -0.075 0.94024
## V76 815.287 8376.102 0.097 0.92286
## V77 1081.988 6709.441 0.161 0.87255
## V78 -1732.606 7698.691 -0.225 0.82287
## V79 475.577 11584.358 0.041 0.96742
## V80 -17252.630 11692.665 -1.476 0.14647
## V81 22885.356 15421.624 1.484 0.14422
## V82 -3424.020 19094.634 -0.179 0.85843
## V83 -8859.197 20140.797 -0.440 0.66197
## V84 1372.520 21418.345 0.064 0.94917
## V85 19799.461 21577.771 0.918 0.36333
## V86 -8217.575 23334.948 -0.352 0.72623
## V87 -5918.362 24626.244 -0.240 0.81108
## V88 -12921.218 23123.750 -0.559 0.57885
## V89 25009.181 19748.368 1.266 0.21136
## V90 -20384.934 20112.386 -1.014 0.31578
## V91 7568.018 23472.335 0.322 0.74850
## V92 5493.059 21833.632 0.252 0.80241
## V93 6426.850 16304.674 0.394 0.69516
## V94 -23192.420 16251.527 -1.427 0.15990
## V95 28156.066 16476.665 1.709 0.09381 .
## V96 -32592.237 16723.444 -1.949 0.05705 .
## V97 21143.009 13794.458 1.533 0.13178
## V98 -24846.926 11056.575 -2.247 0.02916 *
## V99 28810.467 11390.400 2.529 0.01469 *
## V100 -9239.152 5607.683 -1.648 0.10584
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.205 on 49 degrees of freedom
## Multiple R-squared: 0.9969, Adjusted R-squared: 0.9907
## F-statistic: 159 on 100 and 49 DF, p-value: < 2.2e-16# Coeficientes del modelo
# ==============================================================================
df_coeficientes <- modelo$coefficients %>%
enframe(name = "predictor", value = "coeficiente")
df_coeficientes %>%
filter(predictor != "(Intercept)") %>%
ggplot(aes(x = predictor, y = coeficiente)) +
geom_col() +
labs(title = "Coeficientes del modelo OLS") +
theme_bw() +
theme(axis.text.x = element_text(size = 5, angle = 45))
El valor \(R^2_{ajustado}\) obtenido es muy alto (0.9967) lo que indica que el modelo es capaz de predecir con gran exactitud el contenido en grasa de las observaciones con las que se ha entrenado. El hecho de que el modelo en conjunto sea significativo (p-value: < 2.2e-16), pero que muy pocos de los predictores lo sean a nivel individual, es indicativo de una posible redundancia entre los predictores (colinealidad).
¿Cómo de bueno es el modelo prediciendo nuevas observaciones que no han participado en el ajuste? Al tratarse de un modelo de regresión, se emplea como métrica el Mean Square Error (MSE).
\[MSE = \frac{1}{n}\displaystyle\sum_{i=1}^n ( \hat{y_i} - y_i)^2\]
# Predicciones de entrenamiento
# ==============================================================================
predicciones_train <- predict(modelo, newdata = datos_train)
# MSE de entrenamiento
# ==============================================================================
training_mse <- mean((predicciones_train - datos_train$fat)^2)
paste("Error (mse) de entrenamiento:", training_mse)## [1] "Error (mse) de entrenamiento: 0.474058465440023"# Predicciones de test
# ==============================================================================
predicciones_test <- predict(modelo, newdata = datos_test)
# MSE de test
# ==============================================================================
test_mse_ols <- mean((predicciones_test - datos_test$fat)^2)
paste("Error (mse) de test:", test_mse_ols)## [1] "Error (mse) de test: 19.4475565603178"El modelo tiene un MSE muy bajo (0.4740585) cuando predice las mismas observaciones con las que se ha entrenado, pero mucho más alto (19.4475566) al predecir nuevas observaciones. Esto significa que el modelo no es útil, ya que el objetivo es aplicarlo para predecir el contenido en grasa de futuras muestras de carne. A este problema se le conoce como overfitting. Una de las causas por las que un modelo puede sufrir overfitting es la incorporación de predictores innecesarios, que no aportan información o que la información que aportan es redundante.
Best subset selection
Para este caso, la estrategia best subset selection queda descartada por el elevado número de predictores (más de 10).
Stepwise Selection
La función step() de paquete stats permite aplicar el proceso de stepwise selection (both, backward, forward) y seleccionar el mejor modelo en base al AIC.
# Creación y entrenamiento del modelo
# ==============================================================================
modelo <- step(
object = lm(formula = fat ~ ., data = datos_train),
direction = "backward",
scope = list(upper = ~., lower = ~1),
trace = FALSE
)##
## Call:
## lm(formula = fat ~ V1 + V2 + V4 + V6 + V9 + V10 + V11 + V12 +
## V13 + V14 + V16 + V17 + V20 + V21 + V22 + V23 + V25 + V27 +
## V28 + V29 + V31 + V32 + V33 + V35 + V36 + V38 + V39 + V40 +
## V41 + V42 + V45 + V46 + V47 + V49 + V52 + V53 + V54 + V55 +
## V56 + V57 + V60 + V63 + V64 + V66 + V67 + V68 + V69 + V71 +
## V72 + V73 + V74 + V80 + V81 + V83 + V85 + V86 + V88 + V89 +
## V90 + V91 + V93 + V94 + V95 + V96 + V97 + V98 + V99 + V100,
## data = datos_train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.02884 -0.43840 0.02715 0.44309 2.01406
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.766 1.986 3.911 0.000190 ***
## V1 8622.132 3035.869 2.840 0.005701 **
## V2 -9880.824 4468.124 -2.211 0.029826 *
## V4 6047.843 3229.443 1.873 0.064715 .
## V6 -15534.090 2850.522 -5.450 5.31e-07 ***
## V9 5427.059 3856.813 1.407 0.163213
## V10 9282.040 5728.156 1.620 0.109028
## V11 -14533.216 8971.174 -1.620 0.109121
## V12 36430.000 13909.158 2.619 0.010520 *
## V13 -35780.374 14861.298 -2.408 0.018330 *
## V14 23994.842 9780.338 2.453 0.016298 *
## V16 -29192.229 5327.040 -5.480 4.68e-07 ***
## V17 14512.255 4480.295 3.239 0.001740 **
## V20 22642.789 7273.425 3.113 0.002559 **
## V21 -67346.201 15368.394 -4.382 3.49e-05 ***
## V22 72893.780 16868.129 4.321 4.37e-05 ***
## V23 -29745.794 9814.002 -3.031 0.003272 **
## V25 5795.495 3499.058 1.656 0.101530
## V27 -12668.824 5413.419 -2.340 0.021731 *
## V28 20340.849 7457.076 2.728 0.007817 **
## V29 -13296.720 5880.640 -2.261 0.026435 *
## V31 -12353.452 9261.040 -1.334 0.185970
## V32 33054.311 12957.577 2.551 0.012625 *
## V33 -20180.343 8883.982 -2.272 0.025767 *
## V35 -7709.195 5258.443 -1.466 0.146503
## V36 10799.976 4127.673 2.616 0.010596 *
## V38 -14656.434 5692.879 -2.575 0.011858 *
## V39 31559.297 10365.699 3.045 0.003142 **
## V40 -29490.457 14525.041 -2.030 0.045606 *
## V41 28410.290 14310.749 1.985 0.050502 .
## V42 -25336.836 7261.249 -3.489 0.000786 ***
## V45 21169.196 4237.287 4.996 3.31e-06 ***
## V46 -9270.035 5509.645 -1.683 0.096321 .
## V47 -6379.738 3236.986 -1.971 0.052152 .
## V49 2734.498 1388.699 1.969 0.052360 .
## V52 4609.208 3905.912 1.180 0.241431
## V53 -26603.736 7489.880 -3.552 0.000641 ***
## V54 50749.507 8508.558 5.965 6.15e-08 ***
## V55 -47320.742 7747.739 -6.108 3.34e-08 ***
## V56 24464.530 4382.265 5.583 3.06e-07 ***
## V57 -8541.352 2052.916 -4.161 7.85e-05 ***
## V60 2842.997 942.509 3.016 0.003417 **
## V63 12321.513 2936.408 4.196 6.91e-05 ***
## V64 -19103.066 3874.104 -4.931 4.28e-06 ***
## V66 20216.562 5194.936 3.892 0.000203 ***
## V67 -11008.347 8525.919 -1.291 0.200320
## V68 -20199.410 7928.671 -2.548 0.012736 *
## V69 24156.472 4594.223 5.258 1.16e-06 ***
## V71 -20704.311 5236.990 -3.953 0.000164 ***
## V72 10936.232 5679.234 1.926 0.057657 .
## V73 11794.842 4389.386 2.687 0.008743 **
## V74 -9007.432 2899.331 -3.107 0.002608 **
## V80 -13925.312 3533.758 -3.941 0.000171 ***
## V81 16869.817 4750.420 3.551 0.000642 ***
## V83 -7716.592 3496.795 -2.207 0.030162 *
## V85 19616.068 5957.100 3.293 0.001471 **
## V86 -11382.571 6358.209 -1.790 0.077156 .
## V88 -18398.058 6802.816 -2.704 0.008336 **
## V89 27496.839 9676.023 2.842 0.005674 **
## V90 -18962.971 10361.605 -1.830 0.070911 .
## V91 8074.153 6931.918 1.165 0.247527
## V93 12016.898 6703.589 1.793 0.076769 .
## V94 -24598.144 9483.898 -2.594 0.011266 *
## V95 26086.328 8826.230 2.956 0.004087 **
## V96 -30083.880 7189.024 -4.185 7.20e-05 ***
## V97 19622.489 6712.914 2.923 0.004491 **
## V98 -21503.599 6988.577 -3.077 0.002853 **
## V99 24208.149 7247.231 3.340 0.001267 **
## V100 -7382.510 3295.793 -2.240 0.027832 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9771 on 81 degrees of freedom
## Multiple R-squared: 0.9967, Adjusted R-squared: 0.9939
## F-statistic: 355.4 on 68 and 81 DF, p-value: < 2.2e-16# Coeficientes del modelo
# ==============================================================================
df_coeficientes <- modelo$coefficients %>%
enframe(name = "predictor", value = "coeficiente")
df_coeficientes %>%
filter(predictor != "(Intercept)") %>%
ggplot(aes(x = predictor, y = coeficiente)) +
geom_col() +
labs(title = "Coeficientes del modelo Stepwise") +
theme_bw() +
theme(axis.text.x = element_text(size = 6, angle = 45))
## [1] "Número de predictores incluidos en el modelo: 69"# Predicciones de entrenamiento
# ==============================================================================
predicciones_train <- predict(modelo, newdata = datos_train)
# MSE de entrenamiento
# ==============================================================================
training_mse <- mean((predicciones_train - datos_train$fat)^2)
paste("Error (mse) de entrenamiento:", training_mse)## [1] "Error (mse) de entrenamiento: 0.515555790718401"# Predicciones de test
# ==============================================================================
predicciones_test <- predict(modelo, newdata = datos_test)
# MSE de test
# ==============================================================================
test_mse_step <- mean((predicciones_test - datos_test$fat)^2)
paste("Error (mse) de test:", test_mse_step)## [1] "Error (mse) de test: 29.3485192529738"El proceso de stepwise selection devuelve como mejor modelo el formado por 69 de los 100 predictores disponibles. Al haber eliminado predictores del modelo, el error de entrenamiento siempre aumenta, en este caso de 0.5 a 0.6, pero el error de test se ha reducido a 29.3485193.
Ridge
El paquete glmnet incorpora toda una serie de funcionalidades para entrenar modelos lineales (regresión y clasificación) con regularización Ridge, Lasso y Elastic Net. La función glmnet() empleada para entrenar los modelos no permite utilizar formulas ~, necesita una matriz x con el valor de los predictores y un vector y con la variable respuesta. Estas matrices pueden crearse de forma rápida con la función model.matrix(), identificando los predictores y generando las variables dummy necesarias en caso de que los predictores sean cualitativos.
El objeto devuelto por la función glmnet() contiene toda la información relevante del modelo entrenado. Con las funciones plot(), print(), coef() y predict() se puede extraer su información de forma eficiente. Para conocer en más detalle este potente paquete visitar glmnet.
# Matrices de entrenamiento y test
# ==============================================================================
x_train <- model.matrix(fat~., data = datos_train)[, -1]
y_train <- datos_train$fat
x_test <- model.matrix(fat~., data = datos_test)[, -1]
y_test <- datos_test$fatEl resultado de un ajuste por ridge regression depende del hiperparametro \(\lambda\) que determina el grado de penalización. Mediante el argumento lambda, se puede especificar un valor concreto o una secuencia de valores. Si no se tiene conocimiento previo de qué valor de \(\lambda\) es el adecuado, abarcar el rango \(10^{10}\) a \(10^{-2}\), que va desde un modelo muy restrictivo que no contiene ningún predictor, hasta uno sin penalización equivalente al ajuste por mínimos cuadrados. La función glmnet() estandariza por defecto las variables antes de realizar el ajuste del modelo.
# Creación y entrenamiento del modelo
# ==============================================================================
# Para obtener un ajuste con regularización Ridge se indica argumento alpha=0.
# Si no se especifica valor de lambda, se selecciona un rango automático.
modelo <- glmnet(
x = x_train,
y = y_train,
alpha = 0,
nlambda = 100,
standardize = TRUE
)glmnet() almacena en una matriz el valor de los coeficientes de regresión para cada valor de \(lambda\). Esto permite acceder, mediante la función coef(), a los coeficientes obtenidos para un determinado valor de \(lambda\) (que haya sido incluido en el rango cuando se han generado los modelos). También permite representar la evolución de los coeficientes a medida que se incrementa \(lambda\).
# Evolución de los coeficientes en función de lambda
# ==============================================================================
regularizacion <- modelo$beta %>%
as.matrix() %>%
t() %>%
as_tibble() %>%
mutate(lambda = modelo$lambda)
regularizacion <- regularizacion %>%
pivot_longer(
cols = !lambda,
names_to = "predictor",
values_to = "coeficientes"
)
regularizacion %>%
ggplot(aes(x = lambda, y = coeficientes, color = predictor)) +
geom_line() +
scale_x_log10(
breaks = trans_breaks("log10", function(x) 10^x),
labels = trans_format("log10", math_format(10^.x))
) +
labs(title = "Coeficientes del modelo en función de la regularización") +
theme_bw() +
theme(legend.position = "none")
Puede verse como, a medida que aumenta el valor de lambda, la regularización es mayor y el valor de los coeficientes se reduce.
Para identificar el valor de \(lambda\) que da lugar al mejor modelo, se puede recurrir a validación cruzada con la función cv.glmnet().
# Evolución del error en función de lambda
# ==============================================================================
set.seed(123)
cv_error <- cv.glmnet(
x = x_train,
y = y_train,
alpha = 0,
nfolds = 10,
type.measure = "mse",
standardize = TRUE
)
plot(cv_error)
El gráfico muestra el Mean Square Error obtenido por validación cruzada para cada valor de \(lambda\) junto con la barra de error correspondiente. Entre la información almacenada en el objeto cv.glmnet se encuentra el valor de \(lambda\) con el que se consigue el menor error (lambda.min) y el mayor valor de \(lambda\) con el que se consigue el modelo más sencillo que se aleja menos de 1 desviación estándar del error mínimo encontrado.
# Mejor valor lambda encontrado
# ==============================================================================
paste("Mejor valor de lambda encontrado:", cv_error$lambda.min)## [1] "Mejor valor de lambda encontrado: 0.671944767949877"# Mejor valor lambda encontrado + 1sd
# ==============================================================================
# Mayor valor de lambda con el que el test-error no se aleja más de 1sd del mínimo.
paste("Mejor valor de lambda encontrado + 1 desviación estándar:", cv_error$lambda.1se)## [1] "Mejor valor de lambda encontrado + 1 desviación estándar: 1.06992609161865"Se entrena de nuevo el modelo, esta vez empleando el mayor valor de \(lambda\) cuyo error está a menos de una desviación típica del mínimo encontrado en la validación cruzada.
# Mejor modelo lambda óptimo + 1sd
# ==============================================================================
modelo <- glmnet(
x = x_train,
y = y_train,
alpha = 0,
lambda = cv_error$lambda.1se,
standardize = TRUE
)# Coeficientes del modelo
# ==============================================================================
df_coeficientes <- coef(modelo) %>%
as.matrix() %>%
as_tibble(rownames = "predictor") %>%
rename(coeficiente = s0)
df_coeficientes %>%
filter(predictor != "(Intercept)") %>%
ggplot(aes(x = predictor, y = coeficiente)) +
geom_col() +
labs(title = "Coeficientes del modelo Ridge") +
theme_bw() +
theme(axis.text.x = element_text(size = 6, angle = 45))
# Predicciones de entrenamiento
# ==============================================================================
predicciones_train <- predict(modelo, newx = x_train)
# MSE de entrenamiento
# ==============================================================================
training_mse <- mean((predicciones_train - y_train)^2)
paste("Error (mse) de entrenamiento:", training_mse)## [1] "Error (mse) de entrenamiento: 27.6682804462286"# Predicciones de test
# ==============================================================================
predicciones_test <- predict(modelo, newx = x_test)
# MSE de test
# ==============================================================================
test_mse_ridge <- mean((predicciones_test - y_test)^2)
paste("Error (mse) de test:", test_mse_ridge)## [1] "Error (mse) de test: 32.3563303373879"Lasso
El proceso para realizar un ajuste mediante Lasso y la identificación del mejor valor de \(lambda\) es equivalente al seguido en el caso de Ridge pero indicando en la función glmnet() que alpha=1.
# Matrices de entrenamiento y test
# ==============================================================================
x_train <- model.matrix(fat~., data = datos_train)[, -1]
y_train <- datos_train$fat
x_test <- model.matrix(fat~., data = datos_test)[, -1]
y_test <- datos_test$fat# Creación y entrenamiento del modelo
# ==============================================================================
# Para obtener un ajuste con regularización Lasso se indica argumento alpha=1.
# Si no se especifica valor de lambda, se selecciona un rango automático.
modelo <- glmnet(
x = x_train,
y = y_train,
alpha = 1,
nlambda = 100,
standardize = TRUE
)# Evolución de los coeficientes en función de lambda
# ==============================================================================
regularizacion <- modelo$beta %>%
as.matrix() %>%
t() %>%
as_tibble() %>%
mutate(lambda = modelo$lambda)
regularizacion <- regularizacion %>%
pivot_longer(
cols = !lambda,
names_to = "predictor",
values_to = "coeficientes"
)
regularizacion %>%
ggplot(aes(x = lambda, y = coeficientes, color = predictor)) +
geom_line() +
scale_x_log10(
breaks = trans_breaks("log10", function(x) 10^x),
labels = trans_format("log10", math_format(10^.x))
) +
labs(title = "Coeficientes del modelo en función de la regularización") +
theme_bw() +
theme(legend.position = "none")
Puede verse como, a medida que aumenta el valor de \(lambda\), la regularización es mayor y más predictores quedan excluidos (su coeficiente es 0).
Para identificar el valor de \(lambda\) que da lugar al mejor modelo, se puede recurrir a validación cruzada con la función cv.glmnet().
# Evolución del error en función de lambda
# ==============================================================================
set.seed(123)
cv_error <- cv.glmnet(
x = x_train,
y = y_train,
alpha = 1,
nfolds = 10,
type.measure = "mse",
standardize = TRUE
)
plot(cv_error)
# Mejor valor lambda encontrado
# ==============================================================================
paste("Mejor valor de lambda encontrado:", cv_error$lambda.min)## [1] "Mejor valor de lambda encontrado: 0.0277648410373884"# Mejor valor lambda encontrado + 1sd
# ==============================================================================
# Mayor valor de lambda con el que el test-error no se aleja más de 1sd del mínimo.
paste("Mejor valor de lambda encontrado + 1 desviación estándar:", cv_error$lambda.1se)## [1] "Mejor valor de lambda encontrado + 1 desviación estándar: 0.0485198482093112"Se entrena de nuevo el modelo, esta vez empleando el mayor valor de \(lambda\) cuyo error está a menos de una desviación típica del mínimo encontrado en la validación cruzada.
# Mejor modelo lambda óptimo + 1sd
# ==============================================================================
modelo <- glmnet(
x = x_train,
y = y_train,
alpha = 1,
lambda = cv_error$lambda.1se,
standardize = TRUE
)# Coeficientes del modelo
# ==============================================================================
df_coeficientes <- coef(modelo) %>%
as.matrix() %>%
as_tibble(rownames = "predictor") %>%
rename(coeficiente = s0)
df_coeficientes %>%
filter(predictor != "(Intercept)") %>%
ggplot(aes(x = predictor, y = coeficiente)) +
geom_col() +
labs(title = "Coeficientes del modelo Lasso") +
theme_bw() +
theme(axis.text.x = element_text(size = 6, angle = 45))
De los 100 predictores disponibles, el modelo final solo incluye 8 .
# Predicciones de entrenamiento
# ==============================================================================
predicciones_train <- predict(modelo, newx = x_train)
# MSE de entrenamiento
# ==============================================================================
training_mse <- mean((predicciones_train - y_train)^2)
paste("Error (mse) de entrenamiento:", training_mse)## [1] "Error (mse) de entrenamiento: 11.7606265784461"# Predicciones de test
# ==============================================================================
predicciones_test <- predict(modelo, newx = x_test)
# MSE de test
# ==============================================================================
test_mse_lasso <- mean((predicciones_test - y_test)^2)
paste("Error (mse) de test:", test_mse_lasso)## [1] "Error (mse) de test: 12.8756493409093"Elastic Net
PCR
Los modelos de regresión basados en componentes principales pueden ajustarse mediante la función pcr() del paquete pls.
# Creación y entrenamiento del modelo
# ==============================================================================
set.seed(123)
# Importante estandarizar las variables indicándolo con el argumento scale
# Indicando validation = CV, se emplea 10-fold-cross-validation para
# identificar el número óptimo de componentes.
modelo_pcr <- pcr(fat ~ ., data = datos_train, scale = TRUE, validation = "CV")El summary del modelo pcr devuelve la estimación del RMSEP (raíz cuadrada del MSE) para cada posible número de componentes introducidas en el modelo. También se muestra el % de varianza explicada acumulada por cada número de componentes. Esta formación también puede extraerse con las funciones MSEP() y RMSEP().
# Evolución del error en función del número de componentes
# ==============================================================================
error_cv <- MSEP(modelo_pcr, estimate = "CV")
n_componentes_optimo <- which.min(error_cv$val)
error_cv <- data.frame(
componentes = seq_along(as.vector(error_cv$val)) - 1,
mse = as.vector(error_cv$val)
)
error_cv %>% head()ggplot(data = error_cv, aes(x = componentes, y = mse)) +
geom_line() +
geom_vline(xintercept = n_componentes_optimo, color = 'red', linetype = 'dashed') +
labs(title = "Error en función de las componentes incluidas") +
theme_bw()
# Mejor número de componentes encontrados
# ==============================================================================
paste("Número de componentes óptimo:", n_componentes_optimo)## [1] "Número de componentes óptimo: 34"Una vez identificado el número óptimo de componentes, se reentrena el modelo indicando este valor.
# Predicciones de entrenamiento
# ==============================================================================
predicciones_train <- predict(modelo, newdata = datos_train)
# MSE de entrenamiento
# ==============================================================================
training_mse <- mean((predicciones_train - datos_train$fat)^2)
paste("Error (mse) de entrenamiento:", training_mse)## [1] "Error (mse) de entrenamiento: 13.5355390920122"# Predicciones de test
# ==============================================================================
predicciones_test <- predict(modelo, newdata = datos_test)
# MSE de test
# ==============================================================================
test_mse_pcr <- mean((predicciones_test - datos_test$fat)^2)
paste("Error (mse) de test:", test_mse_pcr)## [1] "Error (mse) de test: 16.8356116591763"PLS
Los modelos de regresión basados en partial least squares pueden ajustarse mediante la función plsr() del paquete pls.
# Creación y entrenamiento del modelo
# ==============================================================================
set.seed(123)
# Importante estandarizar las variables indicándolo con el argumento scale
# Indicando validation = CV, se emplea 10-fold-cross-validation para
# identificar el número óptimo de componentes.
modelo_plsr <- plsr(fat ~ ., data = datos_train, scale = TRUE, validation = "CV")El summary del modelo plsr devuelve la estimación del RMSEP (raíz cuadrada del MSE) para cada posible número de componentes introducidas en el modelo. También se muestra el % de varianza explicada acumulada por cada número de componentes. Esta formación también puede extraerse con las funciones MSEP() y RMSEP().
# Evolución del error en función del número de componentes
# ==============================================================================
error_cv <- MSEP(modelo_plsr, estimate = "CV")
n_componentes_optimo <- which.min(error_cv$val)
error_cv <- data.frame(
componentes = seq_along(as.vector(error_cv$val)) - 1,
mse = as.vector(error_cv$val)
)
error_cv %>% head()ggplot(data = error_cv, aes(x = componentes, y = mse)) +
geom_line() +
geom_vline(xintercept = n_componentes_optimo, color = 'red', linetype = 'dashed') +
labs(title = "Error en función de las componentes incluidas") +
theme_bw()
# Mejor número de componentes encontrados
# ==============================================================================
paste("Número de componentes óptimo:", n_componentes_optimo)## [1] "Número de componentes óptimo: 20"Una vez identificado el número óptimo de componentes, se entrena de nuevo el modelo con el valor encontrado.
# Predicciones de entrenamiento
# ==============================================================================
predicciones_train <- predict(modelo, newdata = datos_train)
# MSE de entrenamiento
# ==============================================================================
training_mse <- mean((predicciones_train - datos_train$fat)^2)
paste("Error (mse) de entrenamiento:", training_mse)## [1] "Error (mse) de entrenamiento: 15.1159637691117"# Predicciones de test
# ==============================================================================
predicciones_test <- predict(modelo, newdata = datos_test)
# MSE de test
# ==============================================================================
test_mse_plsr <- mean((predicciones_test - datos_test$fat)^2)
paste("Error (mse) de test:", test_mse_plsr)## [1] "Error (mse) de test: 18.1979168568982"Comparación
df_comparacion <- data.frame(
modelo = c("ols", "Stepwise", "Ridge", "Lasso", "PCR","PLS"),
mse = c(test_mse_ols, test_mse_step, test_mse_ridge,
test_mse_lasso, test_mse_pcr, test_mse_plsr)
)
ggplot(data = df_comparacion, aes(x = modelo, y = mse)) +
geom_col(width = 0.5) +
geom_text(aes(label = round(mse, 2)), vjust = -0.1) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
En este caso el mejor modelo se obtiene aplicando regularización Lasso. De esta forma, se consigue un modelo con mayor poder predictivo y que requiere de menos predictores.
Anexos
Métricas de validación
Cuando se quiere elegir el mejor de entre un conjunto de modelos, es necesario compararlos empleando una métrica capaz de estimar el comportamiento del modelo cuando predice observaciones que con las que no se ha entrenado (error de test). Existen dos tipos de aproximaciones para estimar este error:
Estimación indirecta a partir de una corrección aplicada al error de entrenamiento que compense el bias por overfitting: \(C_p\), \(AIC\), \(BIC\) y \(R^2_{ajustado}\).
Estimación directa del mediante validación simple o validación cruzada.
Mallow’s (Cp)
El estadístico \(C_p\) introduce una penalización (en sentido creciente) de \(2d\hat{\sigma}^2\) a la suma de residuos cuadrados de entrenamiento (RSS) con el objetivo de compensar el hecho de que, a medida que se introducen más predictores, el error de entrenamiento siempre estima a la baja el error de test. Cuanto menor es el valor de \(C_p\), mejor el modelo.
\[C_p= \frac{1}{n}(RSS+2d\hat{\sigma}^2)\]
siendo \(d\) el número de predictores y \(\hat{\sigma}^2\) la estimación de la varianza del error \(\epsilon\).
Dado que requiere estimar \(\hat{\sigma}^2\), este método no es aplicable si el número de predictores es mayor que el número de observaciones. Tampoco es recomendable si el número de predictores se aproxima al de observaciones.
AIC
El estadístico Akaike Information Criterion (AIC) se puede aplicar a una gran cantidad de modelos ajustados mediante maximum likelihood (mínimos cuadrados es un caso particular de maximum likelihood). Para regresión lineal por mínimos cuadrados, el valor de \(AIC\) es proporcional al de \(C_p\) por lo que ambos llevan a la selección del mismo modelo.
\[AIC= \frac{1}{n\hat{\sigma}^2}(RSS+2d\hat{\sigma}^2)\]
BIC
El método Bayesian Point of View (BIC) para modelos de mínimos cuadrados se define como:
\[C_p= \frac{1}{n}(RSS+log(n)d\hat{\sigma}^2)\]
A diferencia de \(C_p\), en \(BIC\) la penalización está determinada por \(log(n)\), siendo n el número de observaciones. Esto implica que, cuando \(n>7\), el método \(BIC\) introduce mayores penalizaciones, tendiendo a seleccionar modelos con menos predictores que los seleccionados por \(C_p\) (y también que \(AIC\)).
R2-ajustado
El valor de \(R^2\) se define como el porcentaje de varianza de la variable respuesta explicada por el modelo respecto del total de varianza observada. En los modelos múltiples, cuantos más predictores se incluyan en el modelo, mayor es el valor de \(R^2\), ya que, por poco que sea, cada predictor va a explicar una parte de la varianza observada. La idea detrás de \(R^2_{ajustado}\) es que, una vez que los predictores correctos se han incluido en el modelo, la varianza extra que se consigue explicar añadiendo más predictores no compensa la penalización.
\[R^{2}_{ajustado} = 1 - \frac{RSS}{TSS}x\frac{n-1}{n-k-1} = R^{2} - (1-R^{2})\frac{n-1}{n-k-1} = 1-\frac{RSS/df_{e}}{TSS/df_{t}}\]
Al contrario que \(C_p\), \(AIC\) y \(BIC\), cuanto mayor sea el valor de \(R^2_{ajustado}\), mejor el modelo. \(R^2_{ajustado}\) no requiere de la estimación de \(\hat{\sigma}^2\) por lo que se puede emplear cuando el número de predictores supera al número de observaciones. \(R^2_{ajustado}\) no es generalizable a modelos no lineales.
Ninguno de los métodos mencionados anteriormente, a excepción de \(R^2_{ajustado}\), se debe utilizar cuando el número de predictores se aproxima o supera al número de observaciones, ya que requieren de la estimación de \(\hat{\sigma}^2\) y está no es precisa es este tipo de situaciones.
Validación cruzada (cross validation)
Los métodos de validación, también conocidos como resampling, son estrategias que permiten estimar la capacidad predictiva de los modelos cuando se aplican a nuevas observaciones, haciendo uso únicamente de los datos de entrenamiento. La idea en la que se basan todos ellos es la siguiente: el modelo se ajusta empleando un subconjunto de observaciones del conjunto de entrenamiento y se evalúa (calcular una métrica que mida como de bueno es el modelo, por ejemplo, mse) con las observaciones restantes. Este proceso se repite múltiples veces y los resultados se agregan y promedian. Gracias a las repeticiones, se compensan las posibles desviaciones que puedan surgir por el reparto aleatorio de las observaciones. La diferencia entre métodos suele ser la forma en la que se generan los subconjuntos de entrenamiento/validación.
k-Fold-Cross-Validation (CV)
Las observaciones de entrenamiento se reparten en k folds (conjuntos) del mismo tamaño. El modelo se ajusta con todas las observaciones excepto las del primer fold y se evalúa prediciendo las observaciones del fold que ha quedado excluido, obteniendo así la primera métrica. El proceso se repite k veces, excluyendo un fold distinto en cada iteración. Al final, se generan k valores de la métrica, que se agregan (normalmente con la media y la desviación típica) generando la estimación final de validación.
Leave-One-Out Cross-Validation (LOOCV)
LOOCV es un caso especial de k-Fold-Cross-Validation en el que el número k de folds es igual al número de observaciones disponibles en el conjunto de entrenamiento. El modelo se ajusta cada vez con todas las observaciones excepto una, que se emplea para evaluar el modelo. Este método supone un coste computacional muy elevado, el modelo se ajusta tantas veces como observaciones de entrenamiento, por lo que, en la práctica, no suele compensar emplearlo.
Repeated k-Fold-Cross-Validation (repeated CV)
Es exactamente igual al método k-Fold-Cross-Validation pero repitiendo el proceso completo n veces. Por ejemplo, 10-Fold-Cross-Validation con 5 repeticiones implica a un total de 50 iteraciones ajuste-validación, pero no equivale a un 50-Fold-Cross-Validation.
One standard error rule
La estimación de error de test de cada modelo tiene asociado una desviación. La norma de one-standar-error recomienda elegir como mejor modelo aquel que contenga menos predictores y cuya estimación del error no se aleje más de una desviación estándar del modelo con menor error. La idea es que, si varios modelos tienen errores semejantes, es decir, son aproximadamente igual de buenos, se debe escoger el más simple de ellos (principio de parsimonia).
Información sesión
sesion_info <- devtools::session_info()
dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
)Bibliografía
Linear Models with R by Julian J.Faraway
An Introduction to Statistical Learning: with Applications in R (Springer Texts in Statistics)
Regularization and variable selection via the elastic net, Hui Zou and Trevor Hastie, J. R. Statist. Soc.B (2005)¿Cómo citar este documento?
Selección de predictores: subset selection, ridge, lasso y reducción de dimensionalidad por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/31_seleccion_de_predictores_subset_selection_ridge_lasso_dimension_reduction
¿Te ha gustado el artículo? Tu ayuda es importante
Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊
Este material, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.