Los modelos de regresión tienen como objetivo encontrar modelos matemáticos para predecir una variable respuesta a partir de variables predictoras. Estos modelos se construyen con la premisa de minimizar los errores residuales cuadráticos (\(RSS\)), sin embargo los modelos clásicos funcionan bien cuandos se cumplen unas determinadas hipótesis iniciales.

Como hemos visto anteriormente es habitual encontrar situaciones donde al menos una de las hipótesis iniciales no se cumple, por lo que no es ideal aplicar este tipo de técnicas. Como solución a este problema existe los modelos de Regresión Penalizada que consisten en minimizar el error residual cuadrático añadindo una función que penalice los coeficientes del modelo que tomen valores muy altos, esto es, \[RSS + \lambda \cdot \text{pen}(\beta),\] donde

Cuando mayor sea el parámetro de penalización, mayor será la contracción de los coeficientes, reduciendo así su impacto en el modelo, ayudando a la selección de las variables.

Para el estudio de estas técnicas, utilicemos el conjunto \(\texttt{clima}\) del paquete \(\texttt{datos}\).

library(datos)
dat <- clima |> 
  glimpse()
## Rows: 26,115
## Columns: 15
## $ origen           <chr> "EWR", "EWR", "EWR", "EWR", "EWR", "EWR", "EWR", "EWR…
## $ anio             <int> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013,…
## $ mes              <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ dia              <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ hora             <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 14, 15, 16, 17…
## $ temperatura      <dbl> 39.02, 39.02, 39.02, 39.92, 39.02, 37.94, 39.02, 39.9…
## $ punto_rocio      <dbl> 26.06, 26.96, 28.04, 28.04, 28.04, 28.04, 28.04, 28.0…
## $ humedad          <dbl> 59.37, 61.63, 64.43, 62.21, 64.43, 67.21, 64.43, 62.2…
## $ direccion_viento <dbl> 270, 250, 240, 250, 260, 240, 240, 250, 260, 260, 260…
## $ velocidad_viento <dbl> 10.35702, 8.05546, 11.50780, 12.65858, 12.65858, 11.5…
## $ velocidad_rafaga <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ precipitacion    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ presion          <dbl> 1012.0, 1012.3, 1012.5, 1012.2, 1011.9, 1012.4, 1012.…
## $ visibilidad      <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 1…
## $ fecha_hora       <dttm> 2013-01-01 01:00:00, 2013-01-01 02:00:00, 2013-01-01…

Vamos a depurar los datos, centrándonos en un único aeropuerto y eliminando las variables que no tienen relación con el clima, además de eliminar los valores ausentes. En concreto, vamos a considerar las variables que hemos comprobado que sí influyen sobre el punto de rocío.

dat_interes <- dat |> 
  filter(
    origen == "JFK"
  ) |> 
  select(punto_rocio,temperatura,humedad,precipitacion,velocidad_rafaga) |> 
  drop_na() |> 
  glimpse()
## Rows: 1,507
## Columns: 5
## $ punto_rocio      <dbl> 17.96, 17.06, 14.00, 10.94, 8.06, 8.06, 8.06, 8.06, 8…
## $ temperatura      <dbl> 37.94, 37.04, 33.08, 30.02, 26.96, 26.06, 26.06, 24.9…
## $ humedad          <dbl> 44.00, 43.85, 44.92, 44.41, 44.25, 45.93, 45.93, 48.0…
## $ precipitacion    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ velocidad_rafaga <dbl> 24.16638, 25.31716, 24.16638, 29.92028, 35.67418, 29.…

En concreto, vamos a estudiar dos modelos de regresión penalizada:

Ambos modelos se encuentran en el paquete \(\texttt{glmnet}\) utilizando la función \(\texttt{glmnet(alpha,lambda)}\), donde el argumento \(\alpha\) indica el tipo de regresión y el argumento \(\lambda\) el valor del parámetro de penalización.

Para poder aplicar ambas técnicas es necesario construir los conjuntos de entrenamiento y validación, y encontrar el valor de \(\lambda\) óptimo. Para este último se procede a utilizar alguna técnica de validación cruzada, que como hemos visto la óptima es k-Fold CV.

set.seed(1234)
library(rsample)
dat_split <- dat_interes |> 
  initial_split(prop = 0.7)

dat_train <- dat_split |> 
  training() |> 
  glimpse()
## Rows: 1,054
## Columns: 5
## $ punto_rocio      <dbl> 12.02, 51.08, 46.04, 62.96, 8.06, 59.00, 5.00, 48.92,…
## $ temperatura      <dbl> 33.98, 82.94, 69.98, 84.02, 44.06, 60.80, 41.00, 75.9…
## $ humedad          <dbl> 39.72, 33.20, 42.29, 49.22, 22.45, 96.22, 21.97, 38.5…
## $ precipitacion    <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.23, 0.00, 0.00, 0.00,…
## $ velocidad_rafaga <dbl> 24.16638, 17.26170, 23.01560, 23.01560, 29.92028, 23.…
dat_test <- dat_split |> 
  testing() |> 
  glimpse()
## Rows: 453
## Columns: 5
## $ punto_rocio      <dbl> 17.96, 17.06, 14.00, 8.06, 8.06, 8.06, 10.04, 10.04, …
## $ temperatura      <dbl> 37.94, 37.04, 33.08, 26.96, 26.06, 24.08, 30.92, 33.9…
## $ humedad          <dbl> 44.00, 43.85, 44.92, 44.25, 45.93, 49.87, 41.13, 36.3…
## $ precipitacion    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ velocidad_rafaga <dbl> 24.16638, 25.31716, 24.16638, 35.67418, 27.61872, 24.…

Regresión Ridge

Este tipo de regresión penalizada se encuentra en la función \(\texttt{glmnet()}\) cuando el argumento \(\alpha\) toma el valor \(0\). Además, en este modelo los coeficientes de regresión dependen del valor de \(\lambda\) por lo que es necesario estudiar previamente cómo se comportan los coeficientes (o las variables predictoras) en función del valor de \(\lambda\).

Nota: Para poder aplicar la función \(\texttt{glmnet()}\) es necesario introducir la variable respuesta como un vector y las variables predictoras como una matriz.

library(glmnet)
## Cargando paquete requerido: Matrix
## 
## Adjuntando el paquete: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## Loaded glmnet 4.1-8
model_rigde <-  dat_train |> 
  select(-punto_rocio) |> 
  as.matrix() |> 
  glmnet(y = dat_train |> 
           select(punto_rocio) |> pull(),
         alpha = 0,
         lambda = 10^seq(5,-2,length = 50))

model_rigde |> 
  plot(xvar = "lambda")

Como se puede observar en el gráfico, a medida que aumenta el valor de \(\lambda\), los coeficientes de cada variable predictora se contraen a 0. Vamos a utilizar la técnica de k-Fold CV para encontrar el valor de \(\lambda\) que minimiza el \(MSE\), para ello utilizamos la función \(\texttt{cv.glmnet(alpha,nfolds)}\).

CV_rigde <- dat_train |> 
  select(-punto_rocio) |> 
  as.matrix() |> 
  cv.glmnet(y = dat_train |> 
           select(punto_rocio) |> pull(),
         alpha = 0,
         nfolds = 10)

CV_rigde |> 
  plot()

Podemos observar que a medida que aumenta el valor de \(\lambda\), el \(MSE\) también aumenta, por ello vamos a buscar el valor mínimo del \(MSE\) y el valor de \(\lambda\) correspondiente.

Ridge <- tibble(
  MSE.min = min(CV_rigde$cvm),
  RMSE.min = sqrt(MSE.min),
  lambda.min = CV_rigde$lambda.min
) |> 
  glimpse()
## Rows: 1
## Columns: 3
## $ MSE.min    <dbl> 4.748857
## $ RMSE.min   <dbl> 2.179187
## $ lambda.min <dbl> 1.620625

De modo que sabemos que el menor \(MSE\) se encuentra tomando un parámetro de penalización \(\lambda\) = 1.6206.

Ahora sí, podemos aplicar el modelo de Regresión Ridge con este parámetro de penalización.

model_rigde <-  dat_train |> 
  select(-punto_rocio) |> 
  as.matrix() |> 
  glmnet(y = dat_train |> 
           select(punto_rocio) |> pull(),
         alpha = 0,
         lambda = Ridge$lambda.min)

library(broom)
model_rigde |> 
  tidy()
## # A tibble: 5 × 5
##   term              step estimate lambda dev.ratio
##   <chr>            <dbl>    <dbl>  <dbl>     <dbl>
## 1 (Intercept)          1 -31.4      1.62     0.987
## 2 temperatura          1   0.834    1.62     0.987
## 3 humedad              1   0.437    1.62     0.987
## 4 precipitacion        1  -4.78     1.62     0.987
## 5 velocidad_rafaga     1  -0.0629   1.62     0.987

Utilizamos la función \(\texttt{tidy()}\) del paquete \(\texttt{broom}\) para mostrar en formato tabla los coeficientes del modelo. Si comparamos estas estimaciones con las dadas por el modelo de regresión lineal múltiple (mostradas a continuación).

dat_train |> 
  lm(punto_rocio ~ .,
     data = _) |> 
  coef()
##      (Intercept)      temperatura          humedad    precipitacion 
##      -36.9731320        0.9018447        0.4719332      -14.5982974 
## velocidad_rafaga 
##       -0.0434804

Podemos observar que los coeficientes con valores más altos, la regresión Ridge los ha penalizado reduciendo su valor, como por ejemplo el coeficiente asociado a la precipitación. Nota: La interpretación de los coeficientes es equivalente a la dada por los modelos de regresión lineal.

Finalmente, podemos proceder a realizar predicciones con este modelo y determinar las medidas de precisión (\(MSE\) y \(RMSE\)). Las predicciones las haremos con la función \(\texttt{predict.glmnet(model,s,newx)}\) del paquete \(\texttt{glmnet}\), donde \(s\) representa el valor de \(\lambda\) y \(newx\) la matriz de variables predictoras.

dat_predict <- dat_test |> 
  mutate(
    predictions_Ridge = model_rigde |> 
      predict.glmnet(newx = dat_test |> 
                select(-punto_rocio) |> 
                as.matrix(),
              s = Ridge$lambda.min) |> 
      as.numeric()
  ) |> 
  glimpse()
## Rows: 453
## Columns: 6
## $ punto_rocio       <dbl> 17.96, 17.06, 14.00, 8.06, 8.06, 8.06, 10.04, 10.04,…
## $ temperatura       <dbl> 37.94, 37.04, 33.08, 26.96, 26.06, 24.08, 30.92, 33.…
## $ humedad           <dbl> 44.00, 43.85, 44.92, 44.25, 45.93, 49.87, 41.13, 36.…
## $ precipitacion     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ velocidad_rafaga  <dbl> 24.16638, 25.31716, 24.16638, 35.67418, 27.61872, 24…
## $ predictions_Ridge <dbl> 17.911185, 17.022607, 14.259185, 8.137855, 8.627575,…

Respecto a las medidas de precisión:

dat_predict |> 
  summarise(
    MSE = mean((punto_rocio - predictions_Ridge)^2),
    RSME = sqrt(MSE)
  ) |> 
  glimpse()
## Rows: 1
## Columns: 2
## $ MSE  <dbl> 5.161918
## $ RSME <dbl> 2.271986

De manera que las predicciones obtenidas tienen un error medio de \(\pm\) 2.272 ºF respecto a los valores reales.

Finalmente, el principal problema que presenta este tipo de regresión penalizada es que utiliza todos las variables predictoras para el modelo dificultando su interpretación. Por ello, surge la siguiente regresión penalizada.

Regresión Lasso

Este tipo de regresión penalizada se encuentra en la función \(\texttt{glmnet()}\) cuando el argumento \(\alpha\) toma el valor \(1\). Además, cuando \(\lambda\) es suficientemente grande, este modelo fuerza a que algunos de los coeficientes deban ser exactamente igual a 0, permitiendo así una selección de las variables predictoras.

Al igual que en la regresión Ridge, vamos a comprobar cómo se comportan los coeficientes de las variables predictoras en función del valor de \(\lambda\).

library(glmnet)
model_lasso <-  dat_train |> 
  select(-punto_rocio) |> 
  as.matrix() |> 
  glmnet(y = dat_train |> 
           select(punto_rocio) |> pull(),
         alpha = 1,
         lambda = 10^seq(5,-2,length = 50))

model_lasso |> 
  plot(xvar = "lambda")

Como se puede observar en el gráfico, a medida que aumenta el valor de \(\lambda\), los coeficientes de cada variable predictora se contraen a 0. Vamos a utilizar la técnica de k-Fold CV para encontrar el valor de \(\lambda\) que minimiza el \(MSE\), para ello utilizamos la función \(\texttt{cv.glmnet(alpha,nfolds)}\).

CV_lasso <- dat_train |> 
  select(-punto_rocio) |> 
  as.matrix() |> 
  cv.glmnet(y = dat_train |> 
           select(punto_rocio) |> pull(),
         alpha = 1,
         nfolds = 10)

CV_lasso |> 
  plot()

Podemos observar que a medida que aumenta el valor de \(\lambda\), el \(MSE\) también aumenta, por ello vamos a buscar el valor mínimo del \(MSE\) y el valor de \(\lambda\) correspondiente. Además, observamos que también a medida que aumenta el valor del parámetro de penalización se introducen menos variables en el modelo (eje horizontal superior).

Lasso <- tibble(
  MSE.min = min(CV_lasso$cvm),
  RMSE.min = sqrt(MSE.min),
  lambda.min = CV_lasso$lambda.min
) |> 
  glimpse()
## Rows: 1
## Columns: 3
## $ MSE.min    <dbl> 2.911346
## $ RMSE.min   <dbl> 1.706267
## $ lambda.min <dbl> 0.05559505

De modo que sabemos que el menor \(MSE\) se encuentra tomando un parámetro de penalización \(\lambda\) = 0.0556.

Ahora sí, podemos aplicar el modelo de Regresión Lasso con este parámetro de penalización.

model_lasso <-  dat_train |> 
  select(-punto_rocio) |> 
  as.matrix() |> 
  glmnet(y = dat_train |> 
           select(punto_rocio) |> pull(),
         alpha = 1,
         lambda = Lasso$lambda.min)

library(broom)
model_lasso |> 
  tidy()
## # A tibble: 5 × 5
##   term              step estimate lambda dev.ratio
##   <chr>            <dbl>    <dbl>  <dbl>     <dbl>
## 1 (Intercept)          1 -36.9    0.0556     0.992
## 2 temperatura          1   0.900  0.0556     0.992
## 3 humedad              1   0.468  0.0556     0.992
## 4 precipitacion        1 -12.1    0.0556     0.992
## 5 velocidad_rafaga     1  -0.0350 0.0556     0.992

Utilizamos la función \(\texttt{tidy()}\) del paquete \(\texttt{broom}\) para mostrar en formato tabla los coeficientes del modelo. Si comparamos estas estimaciones con las dadas por el modelo de regresión lineal múltiple (mostradas a continuación).

dat_train |> 
  lm(punto_rocio ~ .,
     data = _) |> 
  coef()
##      (Intercept)      temperatura          humedad    precipitacion 
##      -36.9731320        0.9018447        0.4719332      -14.5982974 
## velocidad_rafaga 
##       -0.0434804

Podemos observar que la regresión Lasso ha penalizado minimamente a los coeficientes del modelo reduciendo su valor. Además, ha mantenido el mismo número de variables predictoras, de modo que no se puede eliminar ninguna de las existentes sin empeorar el modelo.

Finalmente, podemos proceder a realizar predicciones con este modelo y determinar las medidas de precisión (\(MSE\) y \(RMSE\)). Las predicciones las haremos con la función \(\texttt{predict.glmnet(model,s,newx)}\) del paquete \(\texttt{glmnet}\), donde \(s\) representa el valor de \(\lambda\) y \(newx\) la matriz de variables predictoras.

dat_predict <- dat_test |> 
  mutate(
    predictions_Lasso = model_lasso |> 
      predict.glmnet(newx = dat_test |> 
                select(-punto_rocio) |> 
                as.matrix(),
              s = Lasso$lambda.min) |> 
      as.numeric()
  ) |> 
  glimpse()
## Rows: 453
## Columns: 6
## $ punto_rocio       <dbl> 17.96, 17.06, 14.00, 8.06, 8.06, 8.06, 10.04, 10.04,…
## $ temperatura       <dbl> 37.94, 37.04, 33.08, 26.96, 26.06, 24.08, 30.92, 33.…
## $ humedad           <dbl> 44.00, 43.85, 44.92, 44.25, 45.93, 49.87, 41.13, 36.…
## $ precipitacion     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ velocidad_rafaga  <dbl> 24.16638, 25.31716, 24.16638, 35.67418, 27.61872, 24…
## $ predictions_Lasso <dbl> 16.967390, 16.047263, 13.025024, 6.802984, 7.060608,…

Respecto a las medidas de precisión:

dat_predict |> 
  summarise(
    MSE = mean((punto_rocio - predictions_Lasso)^2),
    RSME = sqrt(MSE)
  ) |> 
  glimpse()
## Rows: 1
## Columns: 2
## $ MSE  <dbl> 3.353111
## $ RSME <dbl> 1.83115

De manera que las predicciones obtenidas tienen un error medio de \(\pm\) 1.8312 ºF respecto a los valores reales.