Técnicas de Validación

Utilicemos el conjunto \(\texttt{clima}\) del paquete \(\texttt{datos}\) para aplicar las distintas técnicas de validación.

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.

dat_interes <- dat |> 
  filter(
    origen == "JFK"
  ) |> 
  select(-c(origen,anio,mes,dia,hora,direccion_viento,visibilidad,fecha_hora)) |> 
  drop_na() |> 
  glimpse()

## Rows: 1,399
## Columns: 7
## $ temperatura      <dbl> 37.94, 37.04, 33.08, 30.02, 26.96, 26.06, 26.06, 24.9…
## $ punto_rocio      <dbl> 17.96, 17.06, 14.00, 10.94, 8.06, 8.06, 8.06, 8.06, 8…
## $ humedad          <dbl> 44.00, 43.85, 44.92, 44.41, 44.25, 45.93, 45.93, 48.0…
## $ velocidad_viento <dbl> 17.26170, 16.11092, 13.80936, 21.86482, 20.71404, 18.…
## $ velocidad_rafaga <dbl> 24.16638, 25.31716, 24.16638, 29.92028, 35.67418, 29.…
## $ precipitacion    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ presion          <dbl> 1012.1, 1013.2, 1014.3, 1015.3, 1016.5, 1016.1, 1016.…

Bien, ahora vamos a construir un modelo de regresión para poder aplicar las distintas técnicas de validación para determinar las medidas de precisión error cuadrático medio (\(MSE\)) y raíz del error cuadrático medio (\(RMSE\))

En primer lugar, como sabemos que el punto de rocío es la temperatura de condensación del vapor de agua que depende de la temperatura atmosférica, la humedad, precipitación, presión atmosférica y viento, vamos a proponer un modelo de regresión lineal múltiple empleando la metodología forward para la selección de variables predictoras.

model.null <- dat_interes |> 
  lm(punto_rocio ~ 1,
     data = _) 

model.full <- dat_interes |> 
  lm(punto_rocio ~ .,
     data = _)

model <- model.null |> 
  step(direction = "forward",
       scope = list(lower = model.null,
                    upper = model.full))

## Start:  AIC=8210.96
## punto_rocio ~ 1
## 
##                    Df Sum of Sq    RSS    AIC
## + temperatura       1    376956 117516 6202.7
## + humedad           1    178624 315848 7585.9
## + presion           1     39223 455249 8097.3
## + velocidad_rafaga  1     17280 477192 8163.2
## + velocidad_viento  1     12050 482423 8178.4
## + precipitacion     1     11569 482903 8179.8
## <none>                          494472 8211.0
## 
## Step:  AIC=6202.71
## punto_rocio ~ temperatura
## 
##                    Df Sum of Sq    RSS    AIC
## + humedad           1    113431   4085 1505.2
## + precipitacion     1      8662 108854 6097.6
## + presion           1      6228 111288 6128.5
## + velocidad_viento  1      1428 116088 6187.6
## + velocidad_rafaga  1       393 117123 6200.0
## <none>                          117516 6202.7
## 
## Step:  AIC=1505.16
## punto_rocio ~ temperatura + humedad
## 
##                    Df Sum of Sq    RSS    AIC
## + precipitacion     1   297.946 3787.2 1401.2
## + velocidad_rafaga  1    82.860 4002.3 1478.5
## + velocidad_viento  1    63.976 4021.2 1485.1
## + presion           1    40.642 4044.5 1493.2
## <none>                          4085.1 1505.2
## 
## Step:  AIC=1401.21
## punto_rocio ~ temperatura + humedad + precipitacion
## 
##                    Df Sum of Sq    RSS    AIC
## + velocidad_rafaga  1    73.059 3714.1 1376.0
## + velocidad_viento  1    60.167 3727.0 1380.8
## + presion           1    19.617 3767.6 1396.0
## <none>                          3787.2 1401.2
## 
## Step:  AIC=1375.96
## punto_rocio ~ temperatura + humedad + precipitacion + velocidad_rafaga
## 
##                    Df Sum of Sq    RSS    AIC
## <none>                          3714.1 1376.0
## + presion           1   2.91733 3711.2 1376.9
## + velocidad_viento  1   0.46884 3713.7 1377.8

Observamos que el mejor modelo (el que menor criterio de información de Akaike presenta (\(AIC\)) es el modelo \[\texttt{punto_rocio} = \beta_0 + \beta_1 \cdot \texttt{temperatura} + \beta_2 \cdot \texttt{humedad} + \beta_3 \cdot \texttt{precipitacion} + \beta_4 \cdot \texttt{velocidad_rafaga} + \varepsilon\]

summary(model)

## 
## Call:
## lm(formula = punto_rocio ~ temperatura + humedad + precipitacion + 
##     velocidad_rafaga, data = dat_interes)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.1408 -0.7211  0.5790  1.0438  5.2923 
## 
## Coefficients:
##                    Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)      -37.534300   0.286012 -131.233  < 2e-16 ***
## temperatura        0.897404   0.002727  329.113  < 2e-16 ***
## humedad            0.488110   0.002460  198.453  < 2e-16 ***
## precipitacion    -16.945236   1.629441  -10.399  < 2e-16 ***
## velocidad_rafaga  -0.040214   0.007680   -5.236 1.89e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.632 on 1394 degrees of freedom
## Multiple R-squared:  0.9925, Adjusted R-squared:  0.9925 
## F-statistic: 4.605e+04 on 4 and 1394 DF,  p-value: < 2.2e-16

Como podemos observar el contraste global (\(F-statistic\)) proporciona un p-valor significativo con lo cual podemos rechazar la hipótesis nula de que ninguna variable predictora influye sobre la variable respuesta, y observando los contrastes individuales para cada variable predictora, concluimos que todas ellas influyen significativamente (p-valor \(< 0.05\)).

Además, podemos ver que el coeficiente de determinación ajustado (\(\bar{R}^2\)) toma un valor de 0.9925, por lo que el modelo explica el 99.25 % de la variabilidad de la variable respuesta. En cuanto al error residual estándar (\(RSE\)), éste indica que el punto de rocío se desvía, en media, 1.63 ºF respecto a la recta de regresión.

Ahora sí, procedamos a aplicar las distintas técnicas de validación para calcular las medidas de precisión \(MSE\) y \(RMSE\).

Validación Cruzada (Cross-Validation, CV)

Esta técnica consiste en dividir, de manera aleatoria, el conjunto de datos en dos subconjuntos: un conjunto de entrenamiento y un conjunto de validación con una proporción 70/30, respectivamente (o en ocasiones, 80/20). Lo haremos empleando las funciones \(\texttt{initial_split(data,prop)}, \texttt{training()}\) y \(\texttt{testing()}\) del paquete \(\texttt{rsample}\).

set.seed(1234)
library(rsample)

dat_split <- dat_interes |> 
  initial_split(prop = 0.7) 

dat_train <- dat_split |> 
  training() |> 
  glimpse()

## Rows: 979
## Columns: 7
## $ temperatura      <dbl> 33.98, 69.08, 52.52, 69.98, 78.98, 93.92, 55.94, 80.9…
## $ punto_rocio      <dbl> 14.00, 44.06, 25.70, 42.98, 42.08, 62.06, 24.08, 51.9…
## $ humedad          <dbl> 43.33, 40.45, 35.07, 37.64, 26.89, 34.88, 28.93, 36.6…
## $ velocidad_viento <dbl> 24.16638, 10.35702, 21.86482, 17.26170, 18.41248, 13.…
## $ velocidad_rafaga <dbl> 28.76950, 25.31716, 29.92028, 21.86482, 26.46794, 18.…
## $ precipitacion    <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,…
## $ presion          <dbl> 1010.0, 1015.6, 1020.2, 1016.2, 1014.1, 1021.6, 1019.…

dat_test <- dat_split |> 
  testing() |> 
  glimpse()

## Rows: 420
## Columns: 7
## $ temperatura      <dbl> 37.94, 37.04, 33.08, 26.96, 26.06, 24.98, 24.08, 30.9…
## $ punto_rocio      <dbl> 17.96, 17.06, 14.00, 8.06, 8.06, 8.06, 8.06, 10.04, 1…
## $ humedad          <dbl> 44.00, 43.85, 44.92, 44.25, 45.93, 48.03, 49.87, 41.1…
## $ velocidad_viento <dbl> 17.26170, 16.11092, 13.80936, 20.71404, 19.56326, 17.…
## $ velocidad_rafaga <dbl> 24.16638, 25.31716, 24.16638, 35.67418, 27.61872, 21.…
## $ precipitacion    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ presion          <dbl> 1012.1, 1013.2, 1014.3, 1016.5, 1016.4, 1016.5, 1016.…

Con el conjunto de entrenamiento aplicamos el modelo de regresión lineal múltiple encontrado y con el conjunto de validación construiremos predicciones para determinar posteriormente las medidas de precisión.

model <- dat_train |> 
  lm(punto_rocio ~ temperatura + humedad + precipitacion + velocidad_rafaga,
     data = _)

summary(model)

## 
## Call:
## lm(formula = punto_rocio ~ temperatura + humedad + precipitacion + 
##     velocidad_rafaga, data = dat_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.0583 -0.7542  0.6080  1.0798  4.7527 
## 
## Coefficients:
##                    Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)      -37.634067   0.347726 -108.229  < 2e-16 ***
## temperatura        0.897264   0.003344  268.313  < 2e-16 ***
## humedad            0.489864   0.003029  161.705  < 2e-16 ***
## precipitacion    -16.037632   1.917937   -8.362  < 2e-16 ***
## velocidad_rafaga  -0.040385   0.009231   -4.375 1.35e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.689 on 974 degrees of freedom
## Multiple R-squared:  0.9919, Adjusted R-squared:  0.9919 
## F-statistic: 2.992e+04 on 4 and 974 DF,  p-value: < 2.2e-16

Para realizar predicciones utilizamos la función \(\texttt{predict(model,newdata)}\).

dat_predict <- dat_test |> 
  mutate(
    predictions = model |> 
      predict(newdata = dat_test)
  ) |> 
  glimpse()

## Rows: 420
## Columns: 8
## $ temperatura      <dbl> 37.94, 37.04, 33.08, 26.96, 26.06, 24.98, 24.08, 30.9…
## $ punto_rocio      <dbl> 17.96, 17.06, 14.00, 8.06, 8.06, 8.06, 8.06, 10.04, 1…
## $ humedad          <dbl> 44.00, 43.85, 44.92, 44.25, 45.93, 48.03, 49.87, 41.1…
## $ velocidad_viento <dbl> 17.26170, 16.11092, 13.80936, 20.71404, 19.56326, 17.…
## $ velocidad_rafaga <dbl> 24.16638, 25.31716, 24.16638, 35.67418, 27.61872, 21.…
## $ precipitacion    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ presion          <dbl> 1012.1, 1013.2, 1014.3, 1016.5, 1016.4, 1016.5, 1016.…
## $ predictions      <dbl> 16.986161, 16.058670, 13.076135, 6.791927, 7.132683, …

CV <- dat_predict |> 
  summarise(
    MSE = mean((punto_rocio - predictions)^2),
    RMSE = sqrt(MSE)
  ) |> 
  glimpse()

## Rows: 1
## Columns: 2
## $ MSE  <dbl> 2.234244
## $ RMSE <dbl> 1.494739

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

Leave-one-out Cross-Validation (LOOCV)

En esta ocasión, el conjunto de entrenamiento contiene todas las observaciones excepto una que se utiliza como conjunto de validación. Por ello, se repite el proceso para todos los datos excluyendo en cada iteración una observación distinta. Para aplicar esta técnica utilizamos el paquete \(\texttt{boot}\) que permite estimar directamente \(MSE\) utilizando esta técnica a través de la función \(\texttt{cv.glm(data,glmfit=model)}\) construyendo el modelo con la función \(\texttt{glm()}\).

De modo que procedemos a construir el modelo de regresión lineal múltiple y posteriormente determinar las medidas de precisión.

library(boot)
model <- dat_interes |> 
  glm(punto_rocio ~ temperatura + humedad + precipitacion + velocidad_rafaga,
     data = _)

LOOCV_model <- dat_interes |>
  cv.glm(glmfit = model) 

LOOCV <- tibble(
  MSE = LOOCV_model$delta[1],
  RMSE = sqrt(MSE)
  ) |> 
  glimpse()

## Rows: 1
## Columns: 2
## $ MSE  <dbl> 2.697359
## $ RMSE <dbl> 1.642364

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

K-Fold Cross-Validation (K-Fold CV)

En esta ocasión, el conjunto de datos se divide de modo aleatorio en \(k\) grupos de aproximadamente el mismo tamaño, de manera que el conjunto de entrenamiento está formado por \(k-1\) grupos y el conjunto de validación por el grupo restante. Este proceso se repite \(k\) veces empleando un grupo distinto en cada ocasión como conjunto de validación.

Es habitual repetir el proceso \(k = 5,10\) veces. Para aplicar esta técnica se emplean las mismas funciones que en la técnica anterior, con la única diferencia de que a la función \(\texttt{cv.glm()}\) se le debe añadir un nuevo argumento \(K\).

library(boot)
model <- dat_interes |> 
  glm(punto_rocio ~ temperatura + humedad + precipitacion + velocidad_rafaga,
     data = _)

kFCV_model <- dat_interes |>
  cv.glm(glmfit = model,
         K = 10) 

kFCV <- tibble(
  MSE = kFCV_model$delta[1],
  RMSE = sqrt(MSE)
  ) |> 
  glimpse()

## Rows: 1
## Columns: 2
## $ MSE  <dbl> 2.693393
## $ RMSE <dbl> 1.641156

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

Comparación de las Técnicas

Si unificamos el \(RMSE\) obtenido por las tres técnicas de validación se puede observar que la técnica de CV es la que proporciona un menor error, pero todas presentan valores aproximados. Realmente, la técnica de CV no es la óptima, ya que proporciona bastante variabilidad como vamos a mostrar a continuación.

Vamos a proceder a repetir 100 veces cada una de las técnicas y mostrar los Box-plots de los errores \(RSME\) obtenidos en cada una de ellas, junto con el valor medio y su desviación típica.

library(rsample)
iter <- 100
RMSE_CV <- numeric(iter)
t0_CV <- proc.time()
for(i in 1:iter){
  dat_split <- dat_interes |> 
    initial_split(prop = 0.7) 
  
  dat_train <- dat_split |> 
    training()
  
  dat_test <- dat_split |> 
    testing() 
  
  model <- dat_train |> 
    lm(punto_rocio ~ temperatura + humedad + precipitacion + velocidad_rafaga,
       data = _)
  
  dat_predict <- dat_test |> 
    mutate(
      predictions = model |> 
        predict(newdata = dat_test)
    )
  
  RMSE_CV[i] <- dat_predict |> 
    summarise(
      MSE = mean((punto_rocio - predictions)^2),
      RMSE = sqrt(MSE)
    ) |> 
    pull(RMSE)
}
t1_CV <- proc.time()

plot_CV <- tibble(
  RMSE = RMSE_CV
) |> 
  ggplot() + 
  aes(x = 1,
      y = RMSE) +
  geom_violin() +
  geom_jitter(colour = "red") + 
  labs(x = "",
       y = "RMSE",
       title = "CV")

library(boot)
iter <- 100
RMSE_LOOCV <- numeric(iter)
t0_LOOCV <- proc.time()
for(i in 1:iter){

  model <- dat_interes |> 
    glm(punto_rocio ~ temperatura + humedad + precipitacion + velocidad_rafaga,
        data = _)
  
  LOOCV_model <- dat_interes |>
    cv.glm(glmfit = model) 
  
  RMSE_LOOCV[i] <- tibble(
    MSE = LOOCV_model$delta[1],
    RMSE = sqrt(MSE)
  ) |> 
    pull(RMSE)
}
t1_LOOCV <- proc.time()

plot_LOOCV <- tibble(
  RMSE = RMSE_LOOCV
) |> 
  ggplot() + 
  aes(x = 1,
      y = RMSE) +
  geom_violin() +
  geom_jitter(colour = "red") + 
  labs(x = "",
       y = "RMSE",
       title = "LOOCV")

library(boot)
iter <- 100
RMSE_kFCV <- numeric(iter)
t0_kFCV <- proc.time()
for(i in 1:iter){
  model <- dat_interes |> 
    glm(punto_rocio ~ temperatura + humedad + precipitacion + velocidad_rafaga,
        data = _)
  
  kFCV_model <- dat_interes |>
    cv.glm(glmfit = model,
           K = 10) 
  
  RMSE_kFCV[i] <- tibble(
    MSE = kFCV_model$delta[1],
    RMSE = sqrt(MSE)
  ) |> 
    pull(RMSE)
}
t1_kFCV <- proc.time()

plot_kFCV <- tibble(
  RMSE = RMSE_kFCV
) |> 
  ggplot() + 
  aes(x = 1,
      y = RMSE) +
  geom_violin() +
  geom_jitter(colour = "red") + 
  labs(x = "",
       y = "RMSE",
       title = "k-Fold CV",
       subtitle = "k = 10")

Para mostrar los tres gráficos a la vez utilizaremos el paquete \(\texttt{ggpubr}\), en particular la función \(\texttt{ggarange}\) que permite representar gráficos creados con el paquete \(\texttt{ggplot2}\) indicando el número de filas y columnas.

library(ggpubr)
ggarrange(plot_CV,
          plot_LOOCV,
          plot_kFCV,
          nrow = 1,
          ncol = 3)

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

Si nos fijamos en la variabilidad, vemos que efectivamente la técnica LOOCV proporciona la menor variabilidad, y que la técnica k-Fold CV también proporciona menor variabilidad que CV.

tibble(
  CV = sd(RMSE_CV),
  LOOCV = sd(RMSE_LOOCV),
  kFCV = sd(RMSE_kFCV)
) |> 
  glimpse()

## Rows: 1
## Columns: 3
## $ CV    <dbl> 0.06893381
## $ LOOCV <dbl> 7.698323e-16
## $ kFCV  <dbl> 0.004319643

Sin embargo, esta técnica presenta un gran coste computacional, como podemos observar a continuación.

tibble(
  CV = t1_CV - t0_CV,
  LOOCV = t1_LOOCV - t0_LOOCV,
  kFCV = t1_kFCV - t0_kFCV,
)

## # A tibble: 5 × 3
##   CV         LOOCV      kFCV      
##   <proc_tim> <proc_tim> <proc_tim>
## 1 0.816      484.293    3.742     
## 2 0.015       33.561    0.215     
## 3 0.838      522.181    3.977     
## 4 0.000        0.000    0.000     
## 5 0.000        0.000    0.000

Por tanto, en términos generales la técnica de validación k-Fold CV proporciona una menor variabilidad que CV y un menor coste computacional que LOOCV.