Paso 1. Depuración del conjunto de datos
Paso 2. Visualización del conjunto de datos
Paso 3. Aplicación de la Técnica de Regresión: Regresión Lineal Múltiple
Paso 4. Diagnosis del modelo
Paso 5. Medidas de Precisión
- Técnicas de Validación: K-Folds Cross-Validation
Paso 6. Predicción

El objetivo de la Regresión Lineal Múltiple es encontrar un modelo relacional entre una variable respuesta \(Y\) y \(p\) variables predictoras \(X_1,...,X_p\) de forma que sea lineal, esto es, \[Y = \beta_0 + \beta_1 \cdot X_1 + \beta_2 \cdot X_2 + \cdots + \beta_p \cdot X_p + \varepsilon,\] donde \(\beta_0, \beta_1, \beta_2, ..., \beta_p\) son coeficientes constantes (a estimar) que representan la ordenada en el origen (\(\beta_0\)) y los efectos medios (\(\beta_j\)) sobre \(Y\) al aumentar la variable predictora \(X_j\) manteniendo el resto de predictores constantes, y \(\varepsilon\) indica el error aleatorio del modelo, que deberá cumplir unas hipótesis iniciales para que el modelo sea válido.

Vayamos paso a paso para encontrar el modelo lineal óptimo utilizando la técnica de Regresión Lineal Múltiple.

Paso 1. Depuración del conjunto de datos

El primer paso es depurar el conjunto de datos donde se encuentren las variables de interés, entendiendo la depuración como la selección de las variables, tratamiento de datos ausentes, eliminación de outliers, etc.

Centrémonos en 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…

Como podemos observar se dispone de 26.115 observaciones referentes a datos meteorológicos en aeropuertos distintos de New York en el año 2013. Entre las distintas variables se encuentran la temperatura (en grados Fahrenheit) (\(\texttt{temperatura}\)), el punto de rocío (en grados Fahrenheit) (\(\texttt{punto_rocio}\)), la humedad relativa (\(\texttt{humedad}\)) y la presión atmosférica (en milibares) (\(\texttt{presion}\)).

Por tanto, vamos a considerar nuestro conjunto de datos considerando estas cuatro variables respecto a los meses de verano con origen el aeropuerto EWR, e identificando y eliminando los valores ausentes.

dat_interes <- dat |> 
  filter(
    origen == "EWR",
    between(mes,6,8)
  ) |> 
  select(temperatura, punto_rocio, humedad, presion) |> 
  glimpse()

## Rows: 2,201
## Columns: 4
## $ temperatura <dbl> 78.08, 77.00, 75.92, 73.94, 73.04, 71.96, 73.94, 75.92, 78…
## $ punto_rocio <dbl> 64.04, 64.94, 66.02, 66.02, 66.02, 66.02, 66.92, 66.92, 66…
## $ humedad     <dbl> 62.05, 66.36, 71.42, 76.31, 78.65, 81.57, 78.72, 73.68, 66…
## $ presion     <dbl> 1016.4, 1016.2, 1016.0, 1015.8, 1016.1, 1016.0, 1016.2, 10…

dat_interes |> 
  mutate(id_row = row_number()) |> 
  filter(if_any(everything(),is.na))

## # A tibble: 250 × 5
##    temperatura punto_rocio humedad presion id_row
##          <dbl>       <dbl>   <dbl>   <dbl>  <int>
##  1        75.9        71.1    84.9      NA     47
##  2        71.6        68      90.1      NA     48
##  3        72.0        69.1    90.7      NA     49
##  4        71.6        69.8    94.1      NA     55
##  5        71.6        69.8    94.1      NA     56
##  6        71.1        70.0   100        NA     58
##  7        73.9        68      88.4      NA     59
##  8        77          66.2    73.6      NA     61
##  9        77          66.2    69.3      NA     62
## 10        73.4        69.8    88.5      NA     66
## # ℹ 240 more rows

dat_interes <- dat_interes |> 
  drop_na() |> 
  glimpse()

## Rows: 1,951
## Columns: 4
## $ temperatura <dbl> 78.08, 77.00, 75.92, 73.94, 73.04, 71.96, 73.94, 75.92, 78…
## $ punto_rocio <dbl> 64.04, 64.94, 66.02, 66.02, 66.02, 66.02, 66.92, 66.92, 66…
## $ humedad     <dbl> 62.05, 66.36, 71.42, 76.31, 78.65, 81.57, 78.72, 73.68, 66…
## $ presion     <dbl> 1016.4, 1016.2, 1016.0, 1015.8, 1016.1, 1016.0, 1016.2, 10…

Finalmente, debemos elegir cuál es nuestra variable respuesta y cuál nuestras variables predictoras, como sabemos que el punto de rocío es la temperatura en la que se condensa el vapor de agua según la presión atmosférica y la humedad, entonces parece razonable pensar que la variable respuesta \(Y\) es \(\texttt{punto_rocio}\) y las \(p = 3\) variables predictoras son \(\texttt{temperatura, humedad, presion}\), proporcionando el modelo: \[\texttt{punto_rocio} = \beta_0 + \beta_1 \cdot \texttt{temperatura} + \beta_2 \cdot \texttt{humedad} + \beta_3 \cdot \texttt{presion} + \varepsilon\]

Paso 2. Visualización del conjunto de datos

Una vez depurado el conjunto de datos, el siguiente paso es realizar una visualización de los datos para ver si existe dicha relación o no entre la variable respuesta y las variables predictoras de modo independiente. En caso afirmativo, también permite observar de qué tipo es, puesto que si no fuese lineal no tendría sentido aplicar la técnica de regresión lineal. Para ello vamos a utilizar un gráfico de dispersión.

La función \(\texttt{ggpairs()}\) del paquete \(\texttt{GGally}\) permite representar la matriz de correlaciones de las variables de interés. En la parte superior se encuentra el coeficiente de correlación de Pearson (junto con una simbología que representa si es significativa la correlación), en la diagonal se muestra la función de densidad de cada variable, y en la parte inferior se muestra el gráfico de dispersión entre los parres de variables.

library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

dat_interes |> 
  ggpairs() +
  theme_minimal()

Como podemos observar la variable respuesta tiene correlación significativa con las variables predictoras, por tanto es razonable incluirlas en el modelo de regresión lineal. En concreto, la temperatura y la humedad presentan correlación positiva con el punto de rocío, mientras que la presión tiene correlación negativa. Por tanto, cuando aumentan la temperatura y la humedad, el punto de rocío también aumenta, sin embargo cuando aumenta la presión, entonces el punto de rocío disminuye.

Paso 3. Aplicación de la Técnica de Regresión: Regresión Lineal Múltiple

Conjuntos de Entrenamiento y Validación

Como hemos comprobado que existe relación entre las variables, el siguiente paso es dividir el conjunto de datos en el Conjunto de Entrenamiento y el Conjunto de Validación con el objetivo de poder validar el modelo de regresión posteriormente. Consideraremos la regla 70/30 (en ocasiones también se suele considerar 80/20), de modo que el conjunto de entrenamiento estará formado por el 70% de las observaciones, mientras que el 30% restante comprenderá el conjunto de validación.

set.seed(1234) #Permite fijar la semilla de valores aleatorios
library(rsample)
dat.split <- dat_interes |> 
  initial_split(prop = 0.7)
  
dat.train <- dat.split |> 
  training() |> 
  glimpse()

## Rows: 1,365
## Columns: 4
## $ temperatura <dbl> 78.98, 75.92, 87.98, 86.00, 66.02, 91.94, 75.92, 78.98, 84…
## $ punto_rocio <dbl> 60.98, 66.02, 75.02, 71.96, 51.98, 75.92, 71.06, 71.96, 71…
## $ humedad     <dbl> 54.08, 71.42, 65.58, 63.01, 60.51, 59.67, 84.92, 79.12, 65…
## $ presion     <dbl> 1012.0, 1011.3, 1009.2, 1009.5, 1016.5, 1013.9, 1011.1, 10…

dat.test <- dat.split |> 
  testing() |> 
  glimpse()

## Rows: 586
## Columns: 4
## $ temperatura <dbl> 78.08, 77.00, 75.92, 73.04, 73.94, 75.92, 78.98, 86.00, 91…
## $ punto_rocio <dbl> 64.04, 64.94, 66.02, 66.02, 66.92, 66.92, 66.92, 68.00, 64…
## $ humedad     <dbl> 62.05, 66.36, 71.42, 78.65, 78.72, 73.68, 66.58, 55.04, 42…
## $ presion     <dbl> 1016.4, 1016.2, 1016.0, 1016.1, 1016.2, 1016.5, 1016.5, 10…

Como consecuencia utilizaremos el conjunto de entrenamiento para entrenar el modelo y el conjunto de validación para validar el modelo y medir su precisión utilizando las Técnicas de Validación.

Aplicación del Modelo de Regresión Lineal Múltiple

Para aplicar el modelo de regresión utilizamos la función \(\texttt{lm()}\).

model <- dat.train |> 
  lm(punto_rocio ~ temperatura + humedad + presion,
     data = _)

model |> 
  summary()

## 
## Call:
## lm(formula = punto_rocio ~ temperatura + humedad + presion, data = dat.train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4081 -0.6204  0.3898  0.9115  1.3003 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -45.347613   5.631052  -8.053 1.75e-15 ***
## temperatura   0.948585   0.004391 216.005  < 2e-16 ***
## humedad       0.466839   0.002119 220.297  < 2e-16 ***
## presion       0.005218   0.005442   0.959    0.338    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.125 on 1361 degrees of freedom
## Multiple R-squared:  0.9802, Adjusted R-squared:  0.9802 
## F-statistic: 2.245e+04 on 3 and 1361 DF,  p-value: < 2.2e-16

Tras aplicar el modelo de regresión lineal, tenemos que interpretar los resultados obtenidos.

Contraste de hipótesis: \(H_0: \beta_1 = \beta_2 = ... = \beta_p = 0\) vs \(H_1:\) algún \(\beta_j \neq 0\)

Es necesario resolver el contraste global donde la hipótesis nula equivale a decir que ninguna de las variables predictoras \(X_1,...,X_p\) tiene influencia sobre la variable respuesta \(Y\). En contrapartida, la hipótesis alternativa significa que alguna variable predictora \(X_j\) sí tiene influencia sobre la variable respuesta. Este contraste se encuentra en la última fila del resumen, donde el p-valor es \(<2.2e-16\), de forma que si fijamos el nivel de significación \(\alpha\) al 5%, entonces obtenemos un resultado significativo y por tanto podemos rechazar la hipótesis nula, y por tanto al menos existe una variable predictora que sí influye sobre la variable respuesta \(\texttt{punto_rocio}\).

Contraste de hipótesis: \(H_0: \beta_j = 0\) vs \(H_1: \beta_j \neq 0, j = 1,...,p\)

Tras resolver el contraste global, es necesario resolver los contrastes para cada variable predictora, en la columna encontramos el p-valor del contraste, concretamente en las fila que represente la variable predictora. En esta ocasión, disponemos de un p-valor \(<2e-16\), que al ser inferior a \(\alpha\), entonces obtenemos un resultado Significativo, y por tanto Rechazamos la hipótesis nula. Con lo cual, podemos decir que las variables predictoras \(\texttt{temperatura, humedad}\) sí influyen sobre la variable respuesta \(\texttt{punto_rocio}\), sin embargo para la variable \(\texttt{presion}\) se tiene un resultado no significativo y por tanto no influye sobre el punto de rocío. De modo que habrá que eliminarla del modelo.

model <- dat.train |> 
  lm(punto_rocio ~ temperatura + humedad,
     data = _)

model |> 
  summary()

## 
## Call:
## lm(formula = punto_rocio ~ temperatura + humedad, data = dat.train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4146 -0.6229  0.3817  0.8887  1.3062 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -39.962575   0.411444  -97.13   <2e-16 ***
## temperatura   0.947854   0.004325  219.18   <2e-16 ***
## humedad       0.466381   0.002065  225.89   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.125 on 1362 degrees of freedom
## Multiple R-squared:  0.9802, Adjusted R-squared:  0.9802 
## F-statistic: 3.368e+04 on 2 and 1362 DF,  p-value: < 2.2e-16

Ahora sí tenemos un modelo donde las variables predictoras sí influyen sobre la variable respuesta.

Interpretación de los coeficientes

Como hemos comprobado que las variables predictoras sí influyen sobre la variable respuesta, entonces debemos proceder a interpretar los coeficientes del modelo. Para ello utilizamos las funciones \(\texttt{coef()}\) y \(\texttt{confint()}\) para obtener las estimaciones y los intervalos de confianza, respectivamente.

model |> 
  coef()

## (Intercept) temperatura     humedad 
## -39.9625753   0.9478539   0.4663815

model |> 
  confint()

##                   2.5 %      97.5 %
## (Intercept) -40.7697081 -39.1554425
## temperatura   0.9393702   0.9563376
## humedad       0.4623312   0.4704317

Podemos observar que las estimaciones para \(\beta_0, \beta_1\) y \(\beta_2\) son, respectivamente, -39.963, 0.948 y 0.466, además tenemos los intervalos de confianza al nivel de confianza del 95%. Esto significa que:

Si no consideramos la temperatura y humedad, el punto de rocío toma un valor medio de -39.963 ºF con un intervalo de confianza, al 95%, de (-40.77,-39.155) ºF.
Cuando la temperatura aumenta 1 ºF, manteniendo constante la humedad, el punto de rocío aumenta, en media, 0.948 ºF con un intervalo de confianza, al 95%, de (0.939,0.956) ºF.
Cuando la humedad relativa aumenta 1 %, manteniendo constante la temperatura, el punto de rocío aumenta, en media, 0.466 ºF con un intervalo de confianza, al 95%, de (0.462,0.47) ºF.

Selección de variables predictoras

Aunque hemos comprobado qué variables predictoras influyen sobre la variable respuesta, podemos preguntarnos si alguna otra variable del conjunto inicial influye también sobre ella. Para resolver esta pregunta vamos a utilizar el Método de Selección de Predictores considerando una de las tres versiones: forward, backward o stepwise. En particular, utilizaremos el método mixto (stepwise) empleando la función \(\texttt{step(model, direction, scope)}\) donde el argumento \(\texttt{direction}\) toma los valores \(\texttt{forward, backward, both}\) y el argumento \(\texttt{scope}\) permite indicar los modelos límites con \(\texttt{lower}\) y \(\texttt{upper}\).

dat.stepwise <- dat |> 
  filter(
    origen == "EWR",
    between(mes,6,8)
  ) |> 
  select(-c(origen,anio,mes,dia,hora,direccion_viento,visibilidad,fecha_hora)) |> 
  drop_na() |> 
  glimpse()

## Rows: 346
## Columns: 7
## $ temperatura      <dbl> 91.04, 89.96, 91.04, 91.94, 89.06, 75.92, 78.98, 80.9…
## $ punto_rocio      <dbl> 64.94, 66.02, 64.04, 62.06, 62.06, 66.92, 66.92, 66.9…
## $ humedad          <dbl> 42.21, 45.34, 40.90, 37.10, 40.60, 73.68, 66.58, 62.4…
## $ velocidad_viento <dbl> 14.96014, 13.80936, 18.41248, 20.71404, 14.96014, 13.…
## $ velocidad_rafaga <dbl> 21.86482, 20.71404, 24.16638, 28.76950, 20.71404, 21.…
## $ precipitacion    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ presion          <dbl> 1014.0, 1013.3, 1012.8, 1012.7, 1012.2, 1011.6, 1011.…

dat.stepwise_split <- dat.stepwise |> 
  initial_split(prop = 0.7) 

dat.stepwise_train <- dat.stepwise_split |> 
  training() |> 
  glimpse()

## Rows: 242
## Columns: 7
## $ temperatura      <dbl> 66.92, 80.96, 86.00, 87.98, 84.02, 71.06, 91.04, 91.0…
## $ punto_rocio      <dbl> 53.96, 64.94, 71.96, 69.98, 62.06, 48.02, 69.98, 69.9…
## $ humedad          <dbl> 63.08, 58.25, 63.01, 55.31, 47.68, 43.92, 50.23, 50.2…
## $ velocidad_viento <dbl> 8.05546, 10.35702, 13.80936, 16.11092, 11.50780, 16.1…
## $ velocidad_rafaga <dbl> 17.26170, 19.56326, 19.56326, 21.86482, 19.56326, 29.…
## $ precipitacion    <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,…
## $ presion          <dbl> 1017.4, 1011.2, 1023.1, 1021.0, 1010.9, 1011.6, 1024.…

dat.stepwise_test <- dat.stepwise_split |> 
  testing() |> 
  glimpse()

## Rows: 104
## Columns: 7
## $ temperatura      <dbl> 91.04, 91.04, 84.92, 87.08, 84.92, 64.94, 59.00, 73.0…
## $ punto_rocio      <dbl> 64.94, 64.04, 68.00, 68.00, 64.04, 42.98, 42.08, 41.0…
## $ humedad          <dbl> 42.21, 40.90, 56.98, 53.18, 49.66, 44.78, 53.36, 31.4…
## $ velocidad_viento <dbl> 14.96014, 18.41248, 13.80936, 8.05546, 17.26170, 12.6…
## $ velocidad_rafaga <dbl> 21.86482, 24.16638, 20.71404, 19.56326, 24.16638, 23.…
## $ precipitacion    <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,…
## $ presion          <dbl> 1014.0, 1012.8, 1010.3, 1010.0, 1008.6, 1013.1, 1015.…

model.null <- dat.stepwise_train |> 
  lm(punto_rocio ~ 1, 
     data = _) #modelo nulo: solo es una constante

model.full <- dat.stepwise_train |> 
  lm(punto_rocio ~ ., 
     data = _) #modelo completo: todas las variables

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

## Start:  AIC=1056.49
## punto_rocio ~ 1
## 
##                    Df Sum of Sq   RSS     AIC
## + temperatura       1    8066.4 10823  923.72
## + humedad           1    2211.4 16678 1028.36
## + presion           1     236.1 18653 1055.45
## <none>                          18889 1056.49
## + velocidad_viento  1     136.9 18752 1056.73
## + precipitacion     1       2.5 18887 1058.46
## + velocidad_rafaga  1       0.7 18889 1058.49
## 
## Step:  AIC=923.72
## punto_rocio ~ temperatura
## 
##                    Df Sum of Sq     RSS     AIC
## + humedad           1   10485.1   337.8   86.70
## + presion           1    1064.9  9758.1  900.65
## + precipitacion     1     300.1 10522.8  918.91
## + velocidad_viento  1     107.0 10716.0  923.31
## <none>                          10822.9  923.72
## + velocidad_rafaga  1      15.9 10807.0  925.36
## - temperatura       1    8066.4 18889.3 1056.49
## 
## Step:  AIC=86.7
## punto_rocio ~ temperatura + humedad
## 
##                    Df Sum of Sq     RSS     AIC
## + precipitacion     1      39.0   298.8   59.03
## <none>                            337.8   86.70
## + velocidad_rafaga  1       1.1   336.7   87.93
## + presion           1       0.8   337.0   88.14
## + velocidad_viento  1       0.2   337.6   88.54
## - humedad           1   10485.1 10822.9  923.72
## - temperatura       1   16340.1 16677.9 1028.36
## 
## Step:  AIC=59.03
## punto_rocio ~ temperatura + humedad + precipitacion
## 
##                    Df Sum of Sq     RSS     AIC
## <none>                            298.8   59.03
## + velocidad_rafaga  1       0.5   298.4   60.66
## + presion           1       0.4   298.5   60.74
## + velocidad_viento  1       0.0   298.8   61.03
## - precipitacion     1      39.0   337.8   86.70
## - humedad           1   10224.0 10522.8  918.91
## - temperatura       1   16117.9 16416.8 1026.54

Podemos observar que el mejor modelo para predecir el punto de rocío es incluyendo, además de las variables predictoras temperatura y humedad, la variable precipitación. Por tanto, este modelo presenta el menor valor de la suma de los cuadrados de los residuos (\(RSS\)) y el menor valor del Criterio de Información de Akaike (\(AIC\)).

Nota: Este criterio nos dice que el modelo que mejor se ajusta a los datos es aquel que tenga un menor valor de \(AIC\).

Con lo cual, vamos a trabajar con el nuevo modelo encontrado: \[\texttt{punto_rocio} = \beta_0 + \beta_1 \cdot \texttt{temperatura} + \beta_2 \cdot \texttt{humedad} + \beta_3 \cdot \texttt{precipitacion} + \varepsilon\]

summary(model.stepwise)

## 
## Call:
## lm(formula = punto_rocio ~ temperatura + humedad + precipitacion, 
##     data = dat.stepwise_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4933 -0.6089  0.3056  0.8127  1.8502 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -44.455773   0.859935 -51.697  < 2e-16 ***
## temperatura     0.967026   0.008535 113.303  < 2e-16 ***
## humedad         0.517823   0.005738  90.240  < 2e-16 ***
## precipitacion -19.012903   3.412691  -5.571 6.82e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.121 on 238 degrees of freedom
## Multiple R-squared:  0.9842, Adjusted R-squared:  0.984 
## F-statistic:  4936 on 3 and 238 DF,  p-value: < 2.2e-16

Vemos que efectivamente, el contraste global es significativo, con lo cual existe al menos una variable predictora que influye sobre la variable respuesta. En concreto, cuando observamos los contrastes individuales, efectivamente todas las variables predictoras influyen significativamente.

Paso 4. Diagnosis del modelo

Tras aplicar e interpretar el modelo es necesario comprobar si es válido, esto es, si se cumplen las hipótesis iniciales del modelo. Estas hipótesis son relativas a los residuos del modelo. Necesitamos recoger los residuos en un objeto tipo \(\texttt{tibble}\), de modo que lo haremos haciendo uso de la función \(\texttt{residuals()}\).

resid <- tibble(resid = residuals(model.stepwise)) |> 
  glimpse()

## Rows: 242
## Columns: 1
## $ resid <dbl> 1.0381603, 0.9422061, 0.6235608, 0.7160860, 0.5764952, 1.0161599…

Relación Lineal entre la variable respuesta y las variables predictoras

Es necesario comprobar si existe una relación lineal entre la variable respuesta y cada una de las variables predictoras. Para ello se representan los residuos en función de cada variable predictora, de manera que se podrá concluir que existe una relación lineal cuando en el gráfico no se observe ningún patrón. En caso de existir algún patrón, una posible solución es utilizar transformaciones a la variable predictora como son \(\log(X), \sqrt{X}, X^2\).

Vamos a añadir los residuos al conjunto de datos de entrenamiento del modelo para aplicar la función \(\texttt{ggpairs()}\), de modo que en su última fila encontremos los gráficos de dispersión entre los residuos y las variables predictoras.

library(ggfortify)
dat.stepwise_train |> 
  select(punto_rocio,temperatura,humedad,precipitacion) |>
  mutate(
    resid = residuals(model.stepwise)
  ) |> 
  ggpairs()

Podemos observar que parece que no existe relación lineal entre los residuos y la variable predictora \(\texttt{humedad}\), pues se ve un patrón claro en forma parabólica. Luego una posible solución puede ser aplicar una transformación.

dat.stepwise_train |> 
  select(humedad) |> 
  mutate(
    humedad.log = log(humedad),
    humedad.sqrt = sqrt(humedad),
    humedad.squared = (humedad)^2,
    resid = residuals(model.stepwise)
  ) |> 
  ggpairs()

Observando la última fila podemos concluir que ninguna transformación elimina el patrón, con lo cual no aplicaremos ninguna transformación a la variable predictora, porque todas ellas implican una correlación con los residuos, algo que no puede ocurrir, pues éstos deben ser independientes a las variables predictoras.

Homocedasticidad de los residuos

Consiste en comprobar si los residuos presentan una varianza constante, esto es, si existe homocedasticidad entre los residuos. Para ello, aplicamos el Test de Breush-Pagan que permite contrastar las hipótesis \[H_0: \texttt{Varianza constante (Homocedasticidad) vs } H_1: \texttt{Varianza no constante (Heterocedasticidad)}\]

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

model.stepwise |> 
  bptest(studentize = FALSE)

## 
##  Breusch-Pagan test
## 
## data:  model.stepwise
## BP = 63.095, df = 3, p-value = 1.281e-13

Como se obtiene un p-valor inferior a 0.05, entonces disponemos de un resultado significativo, y por tanto podemos rechazar la hipótesis nula de homocedasticidad, con lo cual una posible solución será aplicar una transformación a la variable respuesta, como pueden ser: \(\log(Y), \sqrt{Y}\) si la varianza aumenta o \(\exp(Y), Y^2\) si disminuye. Para saber cuál es la que debemos aplicar, vamos a representar el diagnosis.

library(ggfortify)
model.stepwise |> 
  autoplot() +
  theme_minimal()

Aunque el test nos da un resultado significativo, gráficamente se puede observar que no existe ningún patrón para afirmar que la varianza de los residuos no sea constante, aunque es cierto que en el gráfico Scale-Location si parece que existe cierta variabilidad, pero esto puede deberse a la existencia de leverages (que veremos posteriormente).

Independencia de los residuos

Debemos comprobar si los residuos son independientes, para ello se aplica el Test de Durbin-Watson que permite contrastar las hipótesis \[H_0: \texttt{Independencia vs } H_1: \texttt{Dependencia}\]

Este test se encuentra en la función \(\texttt{dwtest(model, alternative = "two.sided")}\) del paquete \(\texttt{lmtest}\).

library(lmtest)
model.stepwise |> 
  dwtest(alternative = "two.sided")

## 
##  Durbin-Watson test
## 
## data:  model.stepwise
## DW = 2.0971, p-value = 0.4467
## alternative hypothesis: true autocorrelation is not 0

Al obtener un p-valor superior a 0.05, entonces no podemos rechazar la hipótesis nula de independencia.

Normalidad de los residuos

Por último también es necesario que los residuos se distribuyan según un modelo normal, por ello se realiza el Test de Kolmogorov-Smirnov-Lilliefors cuyas hipótesis son \[H_0: \texttt{Normalidad vs } H_1: \texttt{No normalidad}\]

Este test se encuentra en la función \(\texttt{lillie.test(residual)}\) del paquete \(\texttt{nortest}\).

library(nortest)
resid |> 
  pull() |> 
  lillie.test()

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  pull(resid)
## D = 0.12595, p-value = 5.647e-10

De nuevo, se obtiene un resultado significativo (p-valor \(<\) 0.05), entonces podemos rechazar la hipótesis nula de normalidad, sin embargo el gráfico Normal Q-Q muestra cómo existen puntos en la parte inferior del gráfico que no se acercan a la diagonal, y de nuevo pueden deberse a la existencia de leverages.

Multicolinealidad

Uno de los principales problemas que puede presentar la regresión múltiple es la existencia de variables predictoras correlacionadas, pues en esa ocasión resulta difícil diferenciar el efecto de cada una de ellas. Para ello es necesario estudiar los factores de inflación de la varianza (\(VIF\)) que indican ausencia de colinealidad cuando toma valores inferiores a 1, pero si toman valores superiores a 5 indican una gran influencia y por tanto un problema. Como solución a esto se propone o bien eliminar dichas variables o bien combinarlas en un único predictor. Para calcular los VIF utilizaremos la función \(\texttt{vif}\) de la librería \(\texttt{car}\).

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

model.stepwise |> 
  vif()

##   temperatura       humedad precipitacion 
##      1.303868      1.367401      1.101965

Como podemos observar todas las variables predictoras presentan valores de \(VIF\) alrededor de 1, con lo cual el modelo no presenta problemas de multicolinealidad. Además, el paquete \(\texttt{performance}\) permite representar los valores de VIF a través de la función \(\texttt{check_collinearity()}\). En el gráfico se puede observar que los valores de inflación de la varianza son bajos.

library(performance)
vif_results <- model.stepwise |> 
  check_collinearity()

vif_results |> 
  plot()

Estudio de valores atípicos

Una cuestión importante en el modelo de regresión es el estudio de valores atípicos para la variable respuesta. Suelen considerarse valores atípicos aquellos que presentan residuos estudentizados fuera del intervalo \([-3,3]\). Vamos a calcular los residuos estudentizados con la función \(\texttt{rstudent(model)}\).

tibble(
  resid.stud = rstudent(model.stepwise)
) |> 
  ggplot() +
  aes(x = seq_along(resid.stud),
      y = abs(resid.stud)) +
  geom_point(aes(color = ifelse(abs(resid.stud) > 3,
                                "red",
                                "black"))) + 
  geom_hline(yintercept = 3,
             color = "red",
             linetype = "dashed") +
  scale_color_identity() +
  labs(x = "",
       y = "|Studentized Residuals|") +
  theme_minimal()

Podemos observar que efectivamente hay valores atípicos en la variable respuesta, de modo que debemos proceder a identificarlos y eliminarlos.

position.outliers <- tibble(
  resid.stud = rstudent(model.stepwise)
) |> 
  mutate(
    id_row = row_number()
  ) |> 
  filter(abs(resid.stud) > 3) |> 
  pull(id_row)

Una vez identificados habrá que eliminarlos de nuestro conjunto de entrenamiento para proceder posteriormente a repetir el modelo.

Estudio de Leverages

Otro de los puntos importantes es la existencia de valores atípicos en las variables predictoras, denonimados leverages. Para determinarlos utilizaremos la función \(\texttt{influence.measures(model)}\) que proporciona una tabla con las observaciones significativamente influyentes en alguna de las variables predictoras del modelo. Se consideran como leverages aquellas observaciones que cumplan que \[hat > 2.5 \cdot (p+1)/n,\] donde \(p\) es el número de variables predictoras y \(n\) el número de observaciones.

infl.leverages <- model.stepwise |> 
  influence.measures() |> 
  summary()

## Potentially influential observations of
##   lm(formula = punto_rocio ~ temperatura + humedad + precipitacion,      data = dat.stepwise_train) :
## 
##     dfb.1_ dfb.tmpr dfb.hmdd dfb.prcp dffit   cov.r   cook.d  hat    
## 12  -0.20   0.71    -1.22_*  18.17_*  18.62_*  7.54_* 76.94_*  0.92_*
## 13   0.01   0.25    -0.68     0.31    -0.97_*  0.80_*  0.22    0.05  
## 48  -0.14   0.26    -0.22     0.14    -0.47_*  0.97    0.06    0.04  
## 55   0.00   0.08    -0.22    -0.13    -0.40_*  0.98    0.04    0.04  
## 65  -0.42   0.30     0.48    -0.07    -0.54_*  0.87_*  0.07    0.03  
## 135 -0.04   0.11    -0.15    -0.21    -0.40_*  1.02    0.04    0.05_*
## 137 -0.12   0.26    -0.31     0.10    -0.58_*  0.93_*  0.08    0.04  
## 148 -0.09   0.23    -0.30    -0.45    -0.83_*  0.93_*  0.17    0.07_*
## 204 -0.28   0.19     0.33    -0.05    -0.37    0.93_*  0.03    0.02  
## 231 -0.08   0.24    -0.37     0.19    -0.61_*  0.92_*  0.09    0.04  
## 233  0.04   0.10    -0.39    -0.08    -0.59_*  0.92_*  0.08    0.04

p <- 3; n <- nrow(dat.stepwise_train)
leverages <- tibble(
  id_row = as.integer(rownames(infl.leverages)),
  hat = infl.leverages[,"hat"]
  ) |>
  filter(hat > 2.5*(p+1)/n) |> 
  print(n = 10)

## # A tibble: 6 × 2
##   id_row    hat
##    <int>  <dbl>
## 1     12 0.918 
## 2     13 0.0492
## 3     48 0.0418
## 4    135 0.0501
## 5    148 0.0679
## 6    231 0.0416

position.leverages <- leverages |> 
  pull(id_row)

De manera que observaciones con alto valor de \(hat\) deben ser eliminadas, si comprobamos cuál es el límite:

hat_limit <- 2.5*(p+1)/n
hat_limit

## [1] 0.04132231

Y volvemos a comprobar la gráfica de diagnosis Residuals vs Leverages vemos que hay observaciones con leverages por encima de 0.0413 que están influyendo bastante en el modelo, y pueden estar provocando esa falta de normalidad y la homocedasticidad que nos indicaban los test de hipótesis.

library(ggfortify)
model.stepwise |> 
  autoplot()

Por todo ello es necesario eliminar los leverages de nuestro conjunto de datos.

Reajuste del modelo

Como hemos observado hay tanto valores atípicos como leverages, por tanto es necesario eliminar dichas observaciones del modelo para poder aplicarlo de nuevo.

position <- union(position.outliers, position.leverages)
dat.stepwise_train <- dat.stepwise_train |> 
  slice(-position)

model.stepwise <- dat.stepwise_train |> 
  lm(punto_rocio ~ temperatura + humedad + precipitacion,
     data = _) 

model.stepwise |> 
  autoplot() +
  theme_minimal()

Al mostrar los gráficos podemos observar que no existe (aparentemente) un patrón claro de variabilidad de los residuos, así como asumir la no normalidad de los residuos. En cuanto a los valores atípicos y leverages.

tibble(
  resid.stud = rstudent(model.stepwise)
) |> 
  ggplot() +
  aes(x = seq_along(resid.stud),
      y = abs(resid.stud)) +
  geom_point(aes(color = ifelse(abs(resid.stud) > 3,
                                "red",
                                "black"))) + 
  geom_hline(yintercept = 3,
             color = "red",
             linetype = "dashed") +
  scale_color_identity() +
  labs(x = "",
       y = "|Studentized Residuals|") +
  theme_minimal()

infl.leverages <- model.stepwise |> 
  influence.measures() |> 
  summary()

## Potentially influential observations of
##   lm(formula = punto_rocio ~ temperatura + humedad + precipitacion,      data = dat.stepwise_train) :
## 
##     dfb.1_ dfb.tmpr dfb.hmdd dfb.prcp dffit   cov.r   cook.d  hat    
## 21  -0.01   0.14    -0.32     0.15    -0.44_*  0.95    0.05    0.03  
## 49   0.19   0.02    -0.63     0.23    -0.72_*  0.86_*  0.12    0.04  
## 52  -0.01   0.08    -0.17     2.97_*   3.08_*  2.42_*  2.32_*  0.62_*
## 78  -0.03   0.21    -0.43     0.21    -0.61_*  0.90_*  0.09    0.04  
## 122 -0.19   0.33    -0.26     0.18    -0.55_*  0.95    0.07    0.04  
## 125 -0.03   0.16    -0.31     0.15    -0.45_*  0.95    0.05    0.03  
## 132 -0.17   0.34    -0.35    -0.46    -0.95_*  0.89_*  0.21    0.07_*
## 174  0.15  -0.04    -0.36     0.13    -0.41_*  0.96    0.04    0.03  
## 198 -0.33   0.21     0.39    -0.06    -0.44_*  0.90_*  0.05    0.02  
## 199  0.16  -0.05    -0.37     0.12    -0.41_*  0.96    0.04    0.03  
## 226  0.02  -0.01    -0.04    -0.73    -0.79_*  1.53_*  0.16    0.35_*

Seguimos teniendo valores atípicos en la variable respuesta y leverages en las variables predictoras, de modo que tendríamos que volver a reiniciar el proceso hasta eliminarlos o éstos presenten valores de hat situados en el límite.

Paso 5. Medidas de Precisión

Tras realizar la diagnosis del modelo comprobando que cumple las hipótesis previas, el siguiente paso es validar el modelo determinando las distintas medidas de precisión como son el Coeficiente de determinación ajustado (\(\bar{R}^2\)) y el Error Estándar Residual (\(RSE\)). Para ello, vamos a utilizar la función \(\texttt{summary(model)}\).

summary(model.stepwise)

## 
## Call:
## lm(formula = punto_rocio ~ temperatura + humedad + precipitacion, 
##     data = dat.stepwise_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0996 -0.5618  0.2760  0.6892  1.3465 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -4.373e+01  7.080e-01 -61.762  < 2e-16 ***
## temperatura    9.495e-01  7.140e-03 132.989  < 2e-16 ***
## humedad        5.333e-01  5.029e-03 106.039  < 2e-16 ***
## precipitacion -1.346e+02  1.903e+01  -7.073 1.79e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9123 on 231 degrees of freedom
## Multiple R-squared:  0.9895, Adjusted R-squared:  0.9894 
## F-statistic:  7257 on 3 and 231 DF,  p-value: < 2.2e-16

En las dos penúltimas líneas del resumen se encuentran las medidas de precisión:

Coeficiente de determinación: \(\texttt{Multiple R-squared}\)
Coeficiente de determinación ajustado: \(\texttt{Adjusted R-squared}\)
Error estándar residual: \(\texttt{Residual standar error}\)

Podemos acceder directamente a ellos:

tibble(
  Coef_Det = summary(model.stepwise)$r.squared,
  Coef_Det_Aj = summary(model.stepwise)$adj.r.squared,
  RSE = summary(model.stepwise)$sigma
) |> 
  print()

## # A tibble: 1 × 3
##   Coef_Det Coef_Det_Aj   RSE
##      <dbl>       <dbl> <dbl>
## 1    0.990       0.989 0.912

De modo que observamos que \(\bar{R}^2\) es superior a 0.7, con lo cual es un modelo óptimo, pues nos indica que el 98.94 % de la variabilidad de la variable \(\texttt{punto_rocio}\) está explicada por el modelo.

Además, \(RSE\) indica que el punto de rocío se desvía, en media, 0.91 ºF respecto a la recta de regresión.

Técnicas de Validación: K-Folds Cross-Validation

Dadas sus ventajas se empleará la validación cruzada con k capas, o lo que es lo mismo, k-Folds Cross-Validation (k-Folds CV). Recordemos que esta técnica divide el conjunto de datos aleatoriamente en k grupos (con aproximadamente tamaños iguales), de forma que utiliza k-1 grupos para el conjunto de entrenamiento y el grupo restante para el conjunto de validación, repitiendo el proceso k veces y eligiendo en cada iteración un grupo distinto como conjunto de validación.

Para ello, emplearemos funciones del paquete \(\texttt{boot}\), en concreto creamos el modelo de regresión lineal con la función \(\texttt{glm()}\) y realizaremos la validación con la función \(\texttt{cv.glm(data,glmfit = model, K)}\) para calcular las medidas de precisión error cuadrático medio (\(MSE\)) y raíz del error cuadrático medio (\(RMSE\)).

library(boot)

## 
## Attaching package: 'boot'

## The following object is masked from 'package:car':
## 
##     logit

dat_interes <- dat |> 
  select(punto_rocio, temperatura, humedad, presion, precipitacion, velocidad_rafaga) |> 
  drop_na()

model.validation <- dat_interes |> 
  glm(punto_rocio ~ .,
      dat = _) 

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

tibble(
  MSE = KFCV$delta[1],
  RMSE = sqrt(MSE)
) |> 
  print()

## # A tibble: 1 × 2
##     MSE  RMSE
##   <dbl> <dbl>
## 1  2.56  1.60

Por tanto, como \(RMSE\) es 1.6, entonces las predicciones tienen un error medio de \(\pm\) 1.6 ºF respecto a los valores reales. En esta ocasión, veamos si este error es pequeño o grande, para ello mostraremos el histograma de la variable respuesta, así como la media y su desviación típica.

dat_interes |> 
  summarise(
    mean = mean(punto_rocio),
    sd = sd(punto_rocio)
  )

## # A tibble: 1 × 2
##    mean    sd
##   <dbl> <dbl>
## 1  30.8  19.7

dat_interes |> 
  ggplot() + 
  aes(x = punto_rocio) +
  geom_histogram(fill = "black",
                 color = "white") + 
  labs(x = "Punto de rocío (ºF)",
       y = "Frecuencia Absoluta",
       title = "Punto de rocío",
       subtitle = "Histograma",
       caption = "Conjunto de datos clima del paquete datos.") +
  theme_minimal()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Tanto su media (y desviación típica) como su histograma muestran que el valor de \(RMSE\) toma un valor alto respecto a los valores de la variable respuesta. Veamos dicha medida normalizada:

dat_interes |> 
  summarize(
    RMSE_norm = sqrt(KFCV$delta[1])/mean(punto_rocio)
  )

## # A tibble: 1 × 1
##   RMSE_norm
##       <dbl>
## 1    0.0519

Al normalizarlo, vemos que \(RSME\) normalizado toma un valor inferior al 10%, entonces no se trata de un error alto y por tanto el modelo es bastante confiable.

Paso 6. Predicción

El objetivo de la regresión lineal es predecir nuevos valores, es por ello que es necesario comprobar las medidas de precisión: error cuadrático medio (\(MSE\)) y raíz del error cuadrático medio (\(RMSE\)) sobre el conjunto de validación.

La predicción de nuevas observaciones se realiza a través de la función \(\texttt{predict(model, newdata)}\).

dat.predict <- dat.stepwise_test |> 
  mutate(
  predictions = predict(model.stepwise, 
                        newdata = dat.stepwise_test)
) |> 
  glimpse()

## Rows: 104
## Columns: 8
## $ temperatura      <dbl> 91.04, 91.04, 84.92, 87.08, 84.92, 64.94, 59.00, 73.0…
## $ punto_rocio      <dbl> 64.94, 64.04, 68.00, 68.00, 64.04, 42.98, 42.08, 41.0…
## $ humedad          <dbl> 42.21, 40.90, 56.98, 53.18, 49.66, 44.78, 53.36, 31.4…
## $ velocidad_viento <dbl> 14.96014, 18.41248, 13.80936, 8.05546, 17.26170, 12.6…
## $ velocidad_rafaga <dbl> 21.86482, 24.16638, 20.71404, 19.56326, 24.16638, 23.…
## $ precipitacion    <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,…
## $ presion          <dbl> 1014.0, 1012.8, 1010.3, 1010.0, 1008.6, 1013.1, 1015.…
## $ predictions      <dbl> 65.22579, 64.52714, 67.29202, 67.31631, 63.38810, 41.…

A continuación, calculamos las medidas de precisión \(MSE\) y \(RMSE\).

dat.predict |> 
  summarise(
    MSE = mean((punto_rocio - predictions)^2),
    RSME = sqrt(MSE)
  ) |> 
  print()

## # A tibble: 1 × 2
##     MSE  RSME
##   <dbl> <dbl>
## 1  6.65  2.58

Tal y como habíamos visto con la técnica de validación, el valor del \(RMSE\) obtenido tras el proceso de predicción, es bajo, de manera que las predicciones obtenidas tienen un error medio de \(\pm\) 2.58 ºF respecto a los valores reales.

Regresión Lineal Múltiple

Mario M. Pizarro

Análisis y Visualización de Datos