Multicolinealidad
La multicolinealidad es un problema que se presenta de forma común en los modelos de regresión lineal en donde se hacer referencia a la asociación lineal entre las variables explicativas del modelo.
Se dice que es un fenómeno muestral porque en la población no se presenta, por ejemplo, se sabe que el consumo puede ser explicado por la riqueza a nivel poblacional, pero cuando se toman muestras ambas variables resultan ser no significativas, ya el valor de sus respectivos estadísticos t son no significativos.
Por lo general se tiende a sospechar de cierto grado de multicolinealidad cuando el valor de coeficiente de determinación R-cuadrado estimado en el modelo es muy alto.
Multicolinealidad perfecta
Suponiendo que se tiene un modelo de regresión ineal con \(k\) variables explicativas \(x_1,x_2,...,x_k\) (donde \(X_1=1\) para todas las observaciones de forma que den cabida al intercepto), se dice que existirá una relación lineal exacta si se satisface la siguiente condición:
\[ \lambda_1x_1+\lambda_2x_2+...+\lambda_kx_k=0 \]
Donde \(\lambda_1,\lambda_2,...\lambda_n\), son constantes que son somultáneamente iguales a cero.
Consecuencias
Los parámetros estimados son indeterminados debido a que no es posible separar las influencias de las distintas variables explicativas debido a que están relacionadas linealmente.
Multicolinealidad Imperfecta
En este caso a la relación lineal entre las variables explicativas se le suma un término denominado error estocástico. Tiene la siguiente definición:
\[ \lambda_1x_1+\lambda_2x_2+...+\lambda_kx_k+v_i=0 \]
Donde \(v_i\) es el término de error estocástico.
Consecuencias
Podemos estimar parámetros por MCO, pero los valores estimados no son muy confiables. Cuando más grande es la correlación, más próximo a cero será la determinante de la matriz \(X´X\) lo cual incrementará las varianzas y covarianzas de los parámetros estimados.
Consecuencias de la Multicolinealidad
Aunque los estimadores MCO son MELI (Mejor Estimador Linealmente Insesgado), presenta varianzas y covarianzas grandes, lo que dificulta su estimación precisa. Debido a la primera consecuencia, los intervalos de confianza tienen a ser mucho más amplios, por lo que se hace propicia una aceptación más fácil de la hipótesis nula (el verdadero coeficiente poblacional). El valor estadístico T de uno o más coeficientes tiende a cero. El coeficiente de determinación R-cuadrado al igual que la prueba F como medidas de asociación conjunta mostrarán valores altos. Los estimadores MCO y sus errores estándar se hacen sensibles a pequeños cambios en los datos.
Detección de la Multicolinealidad
Como la multicolinealidad es en esencia un fenómeno de tipo muestral que surge de información sobre todo no experimental recopilada en la mayoría de las ciencias sociales, no hay un método único para detectarla o medir su fuerza.
R-Cuadrado Alto con Coeficientes no Significativos
Como ya se ha mencionado, es un síntoma “clásico” de multicolinealidad. Si 𝑅2 es alta, es decir, está por encima de 0.8, la prueba 𝐹, en la mayoría de los casos, rechazará la hipótesis de que los coeficientes parciales de pendiente son simultáneamente iguales a cero, pero las pruebas t individuales mostrarán que ningún coeficiente parcial de pendiente, o muy pocos, son estadísticamente diferentes de cero.
Coeficientes de Correlación y Diagramas de Dispersión
Una forma de detectar la alta correlación que puedan existir en cada uno de los predictores es haciendo una matriz de los coeficientes de correlación y observar el valor de cada uno de estos coeficientes. Otra medida para detectar multicolinealidad es hacer dispersiones entre las variables explicativas, y las que tengan una relación lineal muy alta a nivel gráfico, van a generar sospecha de multicolinealidad.
Análisis de Factor inflador de Varianzas (VIF)
Este factor representa la velocidad con la que crecen las varianzas en el modelo de regresión. El VIF muestra la forma como la varianza de un estimador se infla or la presencia de la multicolinealidad. A medida que el grado de colinealidad entre los regresores se acerca a 1, el VIF se acerca a infinito . Es decir, a medida que el grado de colinealidad aumenta, la varianza de un estimador también y, en el límite, se vuelve infinita. Este factor viene dado por la siguiente expresión:
\[ VIF(\beta_i)=\frac{1}{1-R_{xi}^{2}} \]
En el análisis de multicolinealidad se considera que un VIF cercado o mayor a 10 indica la existencia de alta multicolinealidad. Se puede utilizar la siguiente regla para interpretar el factor inflador de varianzas:
| VIF | Estado.de.los.predicotres |
|---|---|
| VIF=1 | No correlacionados |
| 1<VIF>5 | Moderadamente correlacionados |
| VIF>5 a 10 | Altamente correlacionados |
Ejemplo en R
library(haven)
data <- read_dta("C:/Users/daffy/OneDrive/Escritorio/Curso R - ARCHIVOS/1. R BÁSICO/Sesión 8 B - Archivos y tarea/01-ARCHIVOS UTILIZADOS-DATA/DATA-ECONR-BAS-SESION 8-EJEMPLO 3.dta")library(tidyverse)
data<- data %>% select(deudacred,edad,empleo,ingresos)
library(knitr)
library(kableExtra)
data %>% head(15) %>%
kbl(caption = "Datos a utilizar en el modelo",escape = T) %>%
kable_classic_2(html_font = "cambira",full_width=F)| deudacred | edad | empleo | ingresos |
|---|---|---|---|
| 11.359392 | 41 | 17 | 176 |
| 1.362202 | 27 | 10 | 31 |
| 0.856075 | 40 | 15 | 55 |
| 2.658720 | 41 | 15 | 120 |
| 1.787436 | 24 | 2 | 28 |
| 0.392700 | 41 | 5 | 25 |
| 3.833874 | 39 | 20 | 67 |
| 0.128592 | 43 | 12 | 38 |
| 1.358348 | 24 | 3 | 19 |
| 2.777700 | 36 | 0 | 25 |
| 0.182512 | 27 | 0 | 16 |
| 0.252356 | 25 | 4 | 23 |
| 3.929600 | 52 | 24 | 64 |
| 1.715901 | 37 | 6 | 29 |
| 3.703700 | 48 | 22 | 100 |
library(stargazer)
options(scipen = 9999)
modelo_lineal<-lm(deudacred~.,data =data)
stargazer(modelo_lineal,title = "Ejemplo de Regresión Multiple",type = "text")
Ejemplo de Regresión Multiple
===============================================
Dependent variable:
---------------------------
deudacred
-----------------------------------------------
edad 0.0003
(0.009)
empleo 0.019
(0.012)
ingresos 0.028***
(0.002)
Constant 0.084
(0.286)
-----------------------------------------------
Observations 850
R2 0.306
Adjusted R2 0.304
Residual Std. Error 1.774 (df = 846)
F Statistic 124.584*** (df = 3; 846)
===============================================
Note: *p<0.1; **p<0.05; ***p<0.01 edad empleo ingresos
1.503358 1.907924 1.709537 Nota: Los resultados mostrarán que ninguna de las variables especificadas en la regresión generan multicolinealidad en el modelo.
Soluciones al problema de multicolinealidad
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_{ajsutado}^2\)).
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 modelo, esta vez empleando una métrica de validación (validación cruzada,\(C_p\),\(AIC\),\(BIC\) o \(R_{ajsutado}^2\)).
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_{ajsutado}^2\)).
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. Los 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 λ 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\), 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 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(scales)
library(corrr)
# Modelado
# ==============================================================================
library(glmnet)
library(pls)Datos
El set de datos contiene 101 columnas. Las 100 primeras, nombradas como \(V_1,...,V_{100}\) 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).
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.
library(kableExtra)
# Correlación entre columnas numéricas
# ==============================================================================
df_correlaciones <- datos %>%
correlate(method = "pearson") %>%
stretch(remove.dups = TRUE)##
## Correlation method: 'pearson'
## Missing treated using: 'pairwise.complete.obs'df_correlaciones %>% mutate(r_abs = abs(r)) %>% arrange(desc(r_abs)) %>% head(5) %>% kbl() %>% kable_material_dark() %>%
column_spec(column = 3, width = "3cm", background = alpha("blue",0.3))| x | y | r | r_abs |
|---|---|---|---|
| V10 | V11 | 0.9999960 | 0.9999960 |
| V11 | V12 | 0.9999960 | 0.9999960 |
| V9 | V10 | 0.9999958 | 0.9999958 |
| V12 | V13 | 0.9999956 | 0.9999956 |
| V8 | V9 | 0.9999954 | 0.9999954 |
Nota: 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%).
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)
stargazer(modelo,title = "Ejemplo de Regresión Multiple",type = "text")
Ejemplo de Regresión Multiple
===============================================
Dependent variable:
---------------------------
fat
-----------------------------------------------
V1 8,992.864*
(5,135.663)
V2 -10,649.370
(9,342.756)
V3 2,876.109
(14,317.420)
V4 2,367.215
(26,770.080)
V5 1,335.792
(35,389.990)
V6 -16,468.790
(29,401.960)
V7 3,254.574
(17,824.230)
V8 -6,107.216
(11,636.350)
V9 9,222.489
(8,871.189)
V10 10,797.540
(9,583.731)
V11 -17,498.010
(15,182.340)
V12 36,205.860
(27,881.940)
V13 -35,394.630
(37,620.320)
V14 23,245.360
(35,357.100)
V15 8,464.127
(26,419.020)
V16 -37,353.240**
(16,946.590)
V17 13,688.140
(10,577.170)
V18 5,723.458
(10,395.570)
V19 -6,926.451
(15,432.970)
V20 32,790.720
(24,364.410)
V21 -82,403.100**
(33,862.850)
V22 88,160.070**
(36,171.670)
V23 -40,058.490
(26,845.340)
V24 4,262.500
(18,180.180)
V25 5,858.536
(11,999.290)
V26 -188.519
(9,113.771)
V27 -16,091.300
(10,344.800)
V28 25,189.660*
(14,436.960)
V29 -13,930.380
(17,318.150)
V30 -1,236.733
(26,545.480)
V31 -12,272.310
(30,084.360)
V32 29,011.790
(26,980.800)
V33 -11,818.510
(21,170.470)
V34 -5,347.222
(16,151.220)
V35 -10,547.720
(13,533.520)
V36 18,920.500*
(11,164.410)
V37 -8,377.894
(10,157.310)
V38 -8,266.572
(10,876.390)
V39 28,390.690
(17,148.510)
V40 -27,513.870
(26,702.790)
V41 19,956.580
(30,947.110)
V42 -6,916.956
(32,724.990)
V43 -14,498.500
(29,103.310)
V44 -3,703.608
(20,350.040)
V45 34,789.180*
(19,266.420)
V46 -16,613.480
(12,909.350)
V47 -7,777.822
(5,364.359)
V48 2,196.605
(6,629.730)
V49 4,837.740
(10,209.500)
V50 -3,480.071
(14,069.280)
V51 -78.433
(22,816.540)
V52 8,022.376
(26,344.940)
V53 -29,438.110
(23,122.360)
V54 51,184.940***
(17,339.900)
V55 -47,259.800***
(13,783.140)
V56 25,490.800***
(9,004.096)
V57 -8,702.084
(6,652.006)
V58 -1,545.866
(6,229.912)
V59 -234.443
(7,270.375)
V60 3,159.078
(6,288.764)
V61 2,756.859
(5,530.132)
V62 -1,253.152
(6,072.163)
V63 11,834.860
(8,209.067)
V64 -20,446.790*
(11,554.130)
V65 4,869.320
(14,888.730)
V66 14,359.140
(17,656.740)
V67 -8,161.971
(17,957.050)
V68 -20,571.530
(17,588.320)
V69 21,410.850
(14,938.920)
V70 2,546.365
(11,527.680)
V71 -18,824.170
(12,012.390)
V72 5,484.490
(12,190.580)
V73 16,821.860*
(9,147.796)
V74 -10,317.340
(8,521.342)
V75 -657.657
(8,727.037)
V76 815.287
(8,376.102)
V77 1,081.988
(6,709.441)
V78 -1,732.606
(7,698.691)
V79 475.577
(11,584.360)
V80 -17,252.630
(11,692.670)
V81 22,885.360
(15,421.620)
V82 -3,424.020
(19,094.630)
V83 -8,859.197
(20,140.800)
V84 1,372.520
(21,418.350)
V85 19,799.460
(21,577.770)
V86 -8,217.575
(23,334.950)
V87 -5,918.362
(24,626.240)
V88 -12,921.220
(23,123.750)
V89 25,009.180
(19,748.370)
V90 -20,384.930
(20,112.390)
V91 7,568.018
(23,472.330)
V92 5,493.059
(21,833.630)
V93 6,426.850
(16,304.670)
V94 -23,192.420
(16,251.530)
V95 28,156.070*
(16,476.670)
V96 -32,592.240*
(16,723.440)
V97 21,143.010
(13,794.460)
V98 -24,846.930**
(11,056.580)
V99 28,810.470**
(11,390.400)
V100 -9,239.152
(5,607.683)
Constant 7.885**
(3.185)
-----------------------------------------------
Observations 150
R2 0.997
Adjusted R2 0.991
Residual Std. Error 1.205 (df = 49)
F Statistic 159.050*** (df = 100; 49)
===============================================
Note: *p<0.1; **p<0.05; ***p<0.01# 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_minimal() +
theme(axis.text.x = element_text(size = 5, angle = 45))
El valor \(R_{ajustado}^2\) 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()´´´ del 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
)
Ejemplo de Regresión Multiple
===============================================
Dependent variable:
---------------------------
fat
-----------------------------------------------
V1 8,622.132***
(3,035.869)
V2 -9,880.824**
(4,468.124)
V4 6,047.843*
(3,229.443)
V6 -15,534.090***
(2,850.522)
V9 5,427.059
(3,856.813)
V10 9,282.040
(5,728.156)
V11 -14,533.220
(8,971.174)
V12 36,430.000**
(13,909.160)
V13 -35,780.370**
(14,861.300)
V14 23,994.840**
(9,780.338)
V16 -29,192.230***
(5,327.040)
V17 14,512.250***
(4,480.295)
V20 22,642.790***
(7,273.425)
V21 -67,346.200***
(15,368.390)
V22 72,893.780***
(16,868.130)
V23 -29,745.790***
(9,814.002)
V25 5,795.495
(3,499.058)
V27 -12,668.820**
(5,413.419)
V28 20,340.850***
(7,457.076)
V29 -13,296.720**
(5,880.640)
V31 -12,353.450
(9,261.040)
V32 33,054.310**
(12,957.580)
V33 -20,180.340**
(8,883.982)
V35 -7,709.195
(5,258.443)
V36 10,799.980**
(4,127.673)
V38 -14,656.430**
(5,692.879)
V39 31,559.300***
(10,365.700)
V40 -29,490.460**
(14,525.040)
V41 28,410.290*
(14,310.750)
V42 -25,336.840***
(7,261.249)
V45 21,169.200***
(4,237.287)
V46 -9,270.035*
(5,509.645)
V47 -6,379.738*
(3,236.986)
V49 2,734.498*
(1,388.699)
V52 4,609.208
(3,905.912)
V53 -26,603.740***
(7,489.880)
V54 50,749.510***
(8,508.558)
V55 -47,320.740***
(7,747.739)
V56 24,464.530***
(4,382.265)
V57 -8,541.352***
(2,052.916)
V60 2,842.997***
(942.509)
V63 12,321.510***
(2,936.408)
V64 -19,103.070***
(3,874.104)
V66 20,216.560***
(5,194.936)
V67 -11,008.350
(8,525.919)
V68 -20,199.410**
(7,928.671)
V69 24,156.470***
(4,594.223)
V71 -20,704.310***
(5,236.990)
V72 10,936.230*
(5,679.234)
V73 11,794.840***
(4,389.386)
V74 -9,007.432***
(2,899.331)
V80 -13,925.310***
(3,533.758)
V81 16,869.820***
(4,750.420)
V83 -7,716.592**
(3,496.795)
V85 19,616.070***
(5,957.100)
V86 -11,382.570*
(6,358.209)
V88 -18,398.060***
(6,802.816)
V89 27,496.840***
(9,676.023)
V90 -18,962.970*
(10,361.600)
V91 8,074.153
(6,931.918)
V93 12,016.900*
(6,703.589)
V94 -24,598.140**
(9,483.898)
V95 26,086.330***
(8,826.230)
V96 -30,083.880***
(7,189.024)
V97 19,622.490***
(6,712.914)
V98 -21,503.600***
(6,988.577)
V99 24,208.150***
(7,247.231)
V100 -7,382.510**
(3,295.793)
Constant 7.766***
(1.986)
-----------------------------------------------
Observations 150
R2 0.997
Adjusted R2 0.994
Residual Std. Error 0.977 (df = 81)
F Statistic 355.428*** (df = 68; 81)
===============================================
Note: *p<0.1; **p<0.05; ***p<0.01# 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_minimal() +
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.atrix()´´´, 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))
## # A tibble: 8 x 2
## predictor coeficiente
## <chr> <dbl>
## 1 V14 -28.6
## 2 V15 -32.1
## 3 V16 -20.5
## 4 V17 -9.86
## 5 V40 98.4
## 6 V41 34.3
## 7 V51 -31.3
## 8 V52 -17.8De 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"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()## componentes mse
## 1 0 156.42807
## 2 1 123.20740
## 3 2 120.71653
## 4 3 65.24733
## 5 4 19.65839
## 6 5 12.50641ggplot(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()## componentes mse
## 1 0 156.42807
## 2 1 122.67399
## 3 2 69.30160
## 4 3 30.60886
## 5 4 18.29719
## 6 5 10.25280ggplot(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.