Selección de predictores, regularización ridge, lasso, elasticnet y reducción de dimensionalidad
  • Introducción
  • Subset Selection
    • Best subset selection
    • Forward stepwise selection
    • Backward stepwise selection
    • Hibrid (double) Stepwise Selection
  • Regularización (Shrinkage)
    • Ridge
    • Lasso
    • Comparación Ridge y Lasso
    • Elastic net
    • Lambda Search
  • Reducción de dimensionalidad
    • Principal Components Regression PCR
    • Partial least squares (PLS)
  • Ejemplo
    • Paquetes
    • Datos
    • Relación entre variables
    • Modelos
      • Mínimos cuadrados (OLS)
      • Best subset selection
      • Stepwise Selection
      • Ridge
      • Lasso
      • Elastic Net
      • PCR
      • PLS
      • Comparación
  • Anexos
    • Métricas de validación
      • Mallow’s (Cp)
      • AIC
      • BIC
      • R2-ajustado
      • Validación cruzada (cross validation)
      • One standard error rule
  • Información sesión
  • Bibliografía


Más sobre ciencia de datos: cienciadedatos.net

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:

  1. Crear un conjunto de modelos candidatos (todos los posibles o un conjunto considerable de ellos), mediante diferentes combinaciones de los predictores disponibles.

  2. Para cada posible tamaño de modelo (1 predictor, 2 predictores…) se selecciona el mejor basándose en el error de entrenamientor.

  3. Los modelos finalistas de cada tamaño se comparan entre ellos para una métrica de validación (validación cruzada, Cp, AIC, BIC o R2ajustado) 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:

  1. Se genera lo que se conoce como modelo nulo (M0), que es el modelo sin ningún predictor.

  2. 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 (M1).

  3. Se repite el paso anterior para modelos con dos predictores y así sucesivamente hasta llegar al modelo con todos los predictores (Mk).

  4. De entre los mejores modelos seleccionados para cada número de predictores (M0, M1, M2,…,Mk) se identifica el mejor modelo, esta vez empleando una métrica de validación (validación cruzada, Cp, AIC, BIC o R2ajustado). Anexo 1

A pesar de que este método explora todas las posibilidades, tiene dos limitaciones fundamentales:

  • Rrequerimientos computacionales: Se requiere calcular 2p 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 (M0) 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 M1. 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 (M0, M1, M2,…,Mk), se identifica el mejor, esta vez empleando una métrica de validación (validación cruzada, Cp, AIC, BIC o R2ajustado). 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 Mk). 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 (M0). De entre los mejores modelos seleccionados para cada número de predictores (M0, M1, M2,…,Mk) se identifica el mejor, esta vez empleando una métrica de validación (validación cruzada, Cp, AIC, BIC o R2ajustado). 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 (||β||22=∑pj=1β2j). 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 λ. Cuando λ=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 λ aumenta, mayor es la penalización y menor el valor de los predictores.

n∑i=1(yi−β0−p∑j=1βjxij)2+λp∑j=1β2j=suma residuos cuadrados+λp∑j=1β2j

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 λ, 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 λ=∞). 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 (||β||1=∑pj=1|β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 λ. Cuando λ=0, el resultado es equivalente al de un modelo lineal por mínimos cuadrados ordinarios. A medida que λ aumenta, mayor es la penalización y más predictores quedan excluidos.

n∑i=1(yi−β0−p∑j=1βjxij)2+λp∑j=1|βj|=suma residuos cuadrados+λp∑j=1|β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 (αλ||β||1+12(1−α)||β||22). El grado en que influye cada una de las penalizaciones está controlado por el hiperparámetro α. Su valor está comprendido en el intervalo [0,1]. Cuando α=0, se aplica ridge y cuando α=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 (α 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.

∑ni=1(yi−β0−∑pj=1βjxij)22n+λ(αp∑j=1|βj|+1−α2+p∑j=1β2j



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 R2ajustado 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=1nn∑i=1(^yi−yi)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
          )
summary(modelo)
## 
## 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))

paste("Número de predictores incluidos en el modelo:", length(modelo$coefficients))
## [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$fat

El resultado de un ajuste por ridge regression depende del hiperparametro λ 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 λ es el adecuado, abarcar el rango 1010 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))

df_coeficientes %>%
  filter(
    predictor != "(Intercept)",
    coeficiente != 0
)

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.

modelo <- pcr(fat ~ ., data = datos_train, scale = TRUE, ncomp = n_componentes_optimo)
# 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.

modelo <- plsr(fat ~ ., data = datos_train, scale = TRUE, ncomp = n_componentes_optimo)
# 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: Cp, AIC, BIC y R2ajustado.

  • Estimación directa del mediante validación simple o validación cruzada.

Mallow’s (Cp)


El estadístico Cp introduce una penalización (en sentido creciente) de 2dˆσ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 Cp, mejor el modelo.

Cp=1n(RSS+2dˆσ2)

siendo d el número de predictores y ˆσ2 la estimación de la varianza del error ϵ.

Dado que requiere estimar ˆσ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 Cp por lo que ambos llevan a la selección del mismo modelo.

AIC=1nˆσ2(RSS+2dˆσ2)


BIC

El método Bayesian Point of View (BIC) para modelos de mínimos cuadrados se define como:

Cp=1n(RSS+log(n)dˆσ2)

A diferencia de Cp, 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 Cp (y también que AIC).

R2-ajustado

El valor de R2 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 R2, ya que, por poco que sea, cada predictor va a explicar una parte de la varianza observada. La idea detrás de R2ajustado 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.

R2ajustado=1−RSSTSSxn−1n−k−1=R2−(1−R2)n−1n−k−1=1−RSS/dfeTSS/dft

Al contrario que Cp, AIC y BIC, cuanto mayor sea el valor de R2ajustado, mejor el modelo. R2ajustado no requiere de la estimación de ˆσ2 por lo que se puede emplear cuando el número de predictores supera al número de observaciones. R2ajustado no es generalizable a modelos no lineales.

Ninguno de los métodos mencionados anteriormente, a excepción de R2ajustado, 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 ˆσ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 DOI


¿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! 😊


Creative Commons Licence
Este material, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.

Se permite:

  • Compartir: copiar y redistribuir el material en cualquier medio o formato.

  • Adaptar: remezclar, transformar y crear a partir del material.

Bajo los siguientes términos:

  • Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.

  • NoComercial: No puedes utilizar el material para fines comerciales.

  • CompartirIgual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.