Regresión Lineal Simple
Mario M. Pizarro
Análisis y Visualización Datos
El objetivo de la Regresión Lineal Simple es encontrar un modelo relacional entre una variable respuesta \(Y\) y una variable predictora \(X\) de forma que sea lineal, esto es, \[Y = \beta_0 + \beta_1 \cdot X + \varepsilon,\] donde \(\beta_0, \beta_1\) son coeficientes constantes (a estimar) que representan la ordenada en el origen y la pendiente de la recta, respectivamente, 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 Simple.
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{vuelos}\) del paquete \(\texttt{datos}\).
library(datos)
dat <- vuelos |>
glimpse()## Rows: 336,776
## Columns: 19
## $ anio <int> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 201…
## $ mes <int> 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, …
## $ horario_salida <int> 517, 533, 542, 544, 554, 554, 555, 557, 557, 558, 5…
## $ salida_programada <int> 515, 529, 540, 545, 600, 558, 600, 600, 600, 600, 6…
## $ atraso_salida <dbl> 2, 4, 2, -1, -6, -4, -5, -3, -3, -2, -2, -2, -2, -2…
## $ horario_llegada <int> 830, 850, 923, 1004, 812, 740, 913, 709, 838, 753, …
## $ llegada_programada <int> 819, 830, 850, 1022, 837, 728, 854, 723, 846, 745, …
## $ atraso_llegada <dbl> 11, 20, 33, -18, -25, 12, 19, -14, -8, 8, -2, -3, 7…
## $ aerolinea <chr> "UA", "UA", "AA", "B6", "DL", "UA", "B6", "EV", "B6…
## $ vuelo <int> 1545, 1714, 1141, 725, 461, 1696, 507, 5708, 79, 30…
## $ codigo_cola <chr> "N14228", "N24211", "N619AA", "N804JB", "N668DN", "…
## $ origen <chr> "EWR", "LGA", "JFK", "JFK", "LGA", "EWR", "EWR", "L…
## $ destino <chr> "IAH", "IAH", "MIA", "BQN", "ATL", "ORD", "FLL", "I…
## $ tiempo_vuelo <dbl> 227, 227, 160, 183, 116, 150, 158, 53, 140, 138, 14…
## $ distancia <dbl> 1400, 1416, 1089, 1576, 762, 719, 1065, 229, 944, 7…
## $ hora <dbl> 5, 5, 5, 5, 6, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 6, …
## $ minuto <dbl> 15, 29, 40, 45, 0, 58, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5…
## $ fecha_hora <dttm> 2013-01-01 05:00:00, 2013-01-01 05:00:00, 2013-01-…Como podemos observar se dispone de 336.776 observaciones para 19 variables distintas sobre datos de vuelos que despegaron desde distintos aeropuertos de New York en el año 2013. Entre las distintas variables, se encuentran el tiempo de retraso (en minutos) en la salida (\(\texttt{atraso_salida}\)) y el tiempo de retraso (en minutos) en la llegada (\(\texttt{atraso_llegada}\)). En ambas variables, los valores negativos indican salidas o llegadas adelantadas.
Por tanto, vamos a considerar nuestro conjunto de datos considerando solamente los vuelos de la compañía aérea DL que hayan despegado desde el aeropuerto JFK los tres primeros días entre los meses de verano. Reduciremos el conjunto a aquellos vuelos que hayan sufrido retraso en la salida o llegada.
dat_interes <- dat |>
filter(
aerolinea == "DL",
origen == "JFK",
dia <= 3,
between(mes,6,8),
(atraso_salida > 0 | atraso_llegada > 0)
) |>
select(atraso_salida, atraso_llegada) |>
glimpse()## Rows: 338
## Columns: 2
## $ atraso_salida <dbl> 3, 54, 5, 20, -11, 39, 3, 7, 33, 46, 1, 24, 1, 1, 16, 2…
## $ atraso_llegada <dbl> -10, 39, -24, 5, 3, 19, 32, 9, 18, 14, -1, 18, -15, -21…dat_interes |>
mutate(id_row = row_number()) |>
filter(if_any(everything(),is.na))## # A tibble: 0 × 3
## # ℹ 3 variables: atraso_salida <dbl>, atraso_llegada <dbl>, id_row <int>Para determinar qué filas presentan NA vamos a crear una nueva variable que indique el número de fila (\(\texttt{id_row}\)) y posteriormente filtraremos aquellas filas que contengan al menos un NA (\(\texttt{is.na}\)), utilizando la función \(\texttt{if_any()}\) para todas las variables del conjunto de datos \(\texttt{everything()}\). Si existen filas que contienen valores ausentes para alguna de las variables, sería interesante utilizarlas posteriormente en la predicción, pero para el proceso de modelado, vamos a eliminarlas empleando la función \(\texttt{drop_na()}\).
dat_interes <- dat_interes |>
drop_na() |>
glimpse()## Rows: 338
## Columns: 2
## $ atraso_salida <dbl> 3, 54, 5, 20, -11, 39, 3, 7, 33, 46, 1, 24, 1, 1, 16, 2…
## $ atraso_llegada <dbl> -10, 39, -24, 5, 3, 19, 32, 9, 18, 14, -1, 18, -15, -21…Finalmente, debemos elegir cuál es nuestra variable respuesta y cuál nuestra variable predictora, como sabemos que el retraso en la llegada depende del retraso en la salida, entonces la variable respuesta \(Y\) es \(\texttt{atraso_llegda}\) y la variable predictora \(X\) es \(\texttt{atraso_salida}\), esto es, el modelo a buscar tendría la forma: \[\text{atraso_llegada} = \beta_0 + \beta_1 \cdot \text{atraso_salida} + \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, y en caso afirmativo de qué tipo, 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.
dat_interes |>
ggplot() +
aes(x = atraso_salida,
y = atraso_llegada) +
geom_point(col = "blue") +
labs(x = "Tiempo de retraso en la salida (minutos)",
y = "Tiempo de retraso en la llegada (minutos)",
title = "Tiempos de retraso",
subtitle = "Gráfico de dispersión",
caption = "Conjunto de datos vuelos del paquete datos.") +
geom_smooth(method = "lm",
color = "red") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
A simple vista parece posible que existe una relación lineal entre ambas variables, para una mayor coherencia vamos a acompañar dicha afirmación calculando el Coeficiente de Correlación de Pearson entre ambas variables.
dat_interes |>
summarise(
cor.Pearson = cor(atraso_salida,atraso_llegada)
) |>
print()## # A tibble: 1 × 1
## cor.Pearson
## <dbl>
## 1 0.886Podemos observar que el coeficiente de correlación toma un valor superior a 0.7, con lo cual, nos indica que sí existe relación lineal, además esa relación es positiva o directa, esto es, a medida que aumenta una variable, también lo hace la otra. Por tanto, se puede concluir que cuando aumenta la distancia a recorrer, el tiempoo de vuelo también aumenta.
Paso 3. Aplicación de la Técnica de Regresión: Regresión Lineal Simple
Conjuntos de Entrenamiento y Validación
Como hemos comprobado que existe relación entre ambas 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.
Para ello, se utiliza la función \(\texttt{initial_split(data,prop)}\) para generar la información que permita dividir el conjunto, y las funciones \(\texttt{training()}\) y \(\texttt{testing()}\) para construir los conjuntos de entrenamiento y validación, respectivamente, todas ellas se encuentran en el paquete \(\texttt{rsample}\).
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: 236
## Columns: 2
## $ atraso_salida <dbl> 39, 248, 156, 138, 110, 13, 16, -2, -6, 2, 2, 9, 4, 1, …
## $ atraso_llegada <dbl> 24, 231, 136, 98, 102, 6, 22, 15, 24, -14, -29, 26, -23…dat.test <- dat.split |>
testing() |>
glimpse()## Rows: 102
## Columns: 2
## $ atraso_salida <dbl> 3, 54, 5, -11, 7, 1, 2, 10, 0, 1, -5, -2, 99, 33, -3, -…
## $ atraso_llegada <dbl> -10, 39, -24, 3, 9, -15, -17, -8, 36, -15, 31, 25, 84, …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 Simple
Para aplicar el modelo de regresión utilizamos la función \(\texttt{lm()}\).
model <- dat.train |>
lm(atraso_llegada ~ atraso_salida,
data = _)
model |>
summary() ##
## Call:
## lm(formula = atraso_llegada ~ atraso_salida, data = dat.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -51.566 -16.913 0.326 14.690 90.436
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.58819 1.69182 1.53 0.127
## atraso_salida 0.85371 0.02839 30.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.16 on 234 degrees of freedom
## Multiple R-squared: 0.7945, Adjusted R-squared: 0.7936
## F-statistic: 904.6 on 1 and 234 DF, p-value: < 2.2e-16Tras aplicar el modelo de regresión lineal, tenemos que interpretar los resultados obtenidos.
Contraste de hipótesis: \(H_0: \beta_1 = 0\) vs \(H_1: \beta_1 \neq 0\)
En primer lugar, para resolver este contraste debemos fijar el nivel de significación \(\alpha\) al 5%, esto es, \(\alpha = 0.05\). Ahora, en la columna encontramos el p-valor del contraste, concretamente en la 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 la variable predictora \(\texttt{atraso_salida}\) Sí influye sobre la variable respuesta \(\texttt{atraso_llegada}\).
Interpretación de los coeficientes
Como hemos comprobado que la variable predictora sí influye 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) atraso_salida
## 2.5881877 0.8537073model |>
confint()## 2.5 % 97.5 %
## (Intercept) -0.7449519 5.9213272
## atraso_salida 0.7977844 0.9096302Podemos observar que las estimaciones para \(\beta_0\) y \(\beta_1\) son, respectivamente, 2.588 y 0.854, además tenemos los intervalos de confianza al nivel de confianza del 95%. Esto significa que:
Si no consideramos el tiempo de retraso en la salida, el tiempo de retraso en la llegada toma un valor medio de 2.588 minutos con un intervalo de confianza, al 95%, de (-0.745,5.921) minutos.
Cuando el tiempo de retraso en la salida aumenta 1 minuto, el tiempo de retraso en la llegada aumenta, en media, 0.854 minutos con un intervalo de confianza, al 95%, de (0.798,0.91) minutos.
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)) |>
glimpse()## Rows: 236
## Columns: 1
## $ resid <dbl> -11.8827733, 16.6923960, 0.2334698, -22.3997984, 5.5040067, -7.6…Media de los residuos nula
La primera hipótesis es referente a la media de los residuos que es necesario que sea igual a 0.
resid |>
summarise(mean = mean(resid)) |>
print()## # A tibble: 1 × 1
## mean
## <dbl>
## 1 -6.35e-16resid |>
ggplot() +
aes(x = seq_along(resid),
y = resid) +
geom_point(color = "blue") +
geom_hline(yintercept = 0,
color = "red") +
labs(x = "",
y = "Residuos",
title = "Modelo de Regresión Lineal Simple",
subtitle = "Residuos") +
theme_minimal()
Como podemos observar la media de los residuos es muy pequeña, con lo cual se asume que es 0.
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)}\]
En esta ocasión nos interesa un p-valor superior al nivel de significación \(\alpha\), dado que queremos No Rechazar la hipótesis nula, o lo que es lo mismo, obtener un resultado No Significativo.
El Test de Breusch-Pagan se encuentra en la función \(\texttt{bptest(model)}\) del paquete \(\texttt{lmtest}\). Además, añadimos el argumento \(\texttt{studentize} = \texttt{FALSE}\) cuando los residuos no se encuentren en el intervalo [-1,1].
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericmodel |>
bptest(studentize = FALSE)##
## Breusch-Pagan test
##
## data: model
## BP = 2.8017, df = 1, p-value = 0.09417Como se obtiene un p-valor superior a 0.05, entonces disponemos de un resultado no significativo, y por tanto No podemos rechazar la hipótesis nula de homocedasticidad.
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 |>
dwtest(alternative = "two.sided")##
## Durbin-Watson test
##
## data: model
## DW = 2.202, p-value = 0.1203
## alternative hypothesis: true autocorrelation is not 0Al 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.041954, p-value = 0.396De nuevo, como se obtiene un resultado no significativo (p-valor \(>\) 0.05), entonces no podemos rechazar la hipótesis nula de normalidad.
Diagnosis Gráfica de los residuos
Finalmente, vamos a utilizar la función \(\texttt{autoplot()}\) del paquete \(\texttt{ggfortify}\) para representar distintos gráficos para el diagnosis del modelo.
library(ggfortify)
model |>
autoplot() +
theme_minimal()
En este gráfico se divide en 4 paneles:
Residuals vs Fitted: muestra los valores ajustados por el modelo frente a los residuos, este diagrama permite comprobar si existe o no homocedasticidad. Si se observa un patrón donde a medida que aumentan los valores ajustados, los residuos van dispersándose, entonces no presentan homocedasticidad, o equivalentemente, la varianza no es constante.
Scale - Location: muestra los valores ajustados por el modelo frente a la raíz cuadrada de los residuos estandarizados. Mismo argumento que el gráfico anterior.
Normal Q-Q: muestra los cuantiles teóricos de la distribución normal frente a los residuos estandarizados. Se asume que los residuos siguen un modelo normal cuando se encuentran alrededor de la diagonal.
Residuals vs Leverage: muestra los valores que pueden influir en el modelo (leverage) frente a los residuos estandarizados. Los residuos con grandes valores de Leverage implica que dichas observaciones influyen demasiado en el modelo y habría que tratarlos.
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 (\(R^2\)) y el Error Estándar Residual (\(RSE\)). Para ello, vamos a utilizar la función \(\texttt{summary(model)}\).
summary(model)##
## Call:
## lm(formula = atraso_llegada ~ atraso_salida, data = dat.train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -51.566 -16.913 0.326 14.690 90.436
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.58819 1.69182 1.53 0.127
## atraso_salida 0.85371 0.02839 30.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 22.16 on 234 degrees of freedom
## Multiple R-squared: 0.7945, Adjusted R-squared: 0.7936
## F-statistic: 904.6 on 1 and 234 DF, p-value: < 2.2e-16En las dos penúltimas líneas del resumen se encuentran las medidas de precisió:
- Coeficiente de determinación: \(\texttt{Multiple R-squared}\)
- Error estándar residual: \(\texttt{Residual standar error}\)
Podemos acceder directamente a ellos:
tibble(
Coef_Det = summary(model)$r.squared,
RSE = summary(model)$sigma
) |>
print()## # A tibble: 1 × 2
## Coef_Det RSE
## <dbl> <dbl>
## 1 0.794 22.2De modo que observamos que \(R^2\) es superior a 0.7, con lo cual es un modelo óptimo, pues nos indica que el 79.45 % de la variabilidad de la variable \(\texttt{atraso_llegada}\) está explicada por el modelo.
Además, \(RSE\) indica que el tiempo de retraso en la llegada se desvía, en media, 22.16 minutos 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)
model.validation <- dat_interes |>
glm(atraso_llegada ~ atraso_salida,
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 455. 21.3Por tanto, como \(RMSE\) es 21.32, entonces las predicciones tienen un error medio de \(\pm\) 21.32 minutos 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(atraso_llegada),
sd = sd(atraso_llegada)
) ## # A tibble: 1 × 2
## mean sd
## <dbl> <dbl>
## 1 28.0 45.9dat_interes |>
ggplot() +
aes(x = atraso_llegada) +
geom_histogram(fill = "black",
color = "white") +
labs(x = "Tiempo de retraso en la llegada (minutos)",
y = "Frecuencia Absoluta",
title = "Tiempo de retraso en la llegada",
subtitle = "Histograma",
caption = "Conjunto de datos vuelos 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 muy alto respecto a los valores de la variable respuesta. Veamos dicha medida normalizada:
dat_interes |>
summarize(
RMSE_norm = sqrt(KFCV$delta[1])/mean(atraso_llegada)
)## # A tibble: 1 × 1
## RMSE_norm
## <dbl>
## 1 0.762Al normalizarlo, vemos que \(RSME\) normalizado toma un valor superior al 20%, entonces sí que se trata de un error alto y por tanto el modelo es poco confiable, a pesar de tener un coeficiente de determinación alto.
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 <- tibble(
atraso_salida = dat.test$atraso_salida,
atraso_llegada = dat.test$atraso_llegada,
predictions = predict(model,
newdata = dat.test)
) |>
glimpse()## Rows: 102
## Columns: 3
## $ atraso_salida <dbl> 3, 54, 5, -11, 7, 1, 2, 10, 0, 1, -5, -2, 99, 33, -3, -…
## $ atraso_llegada <dbl> -10, 39, -24, 3, 9, -15, -17, -8, 36, -15, 31, 25, 84, …
## $ predictions <dbl> 5.1493096, 48.6883832, 6.8567243, -6.8025929, 8.5641389…A continuación, calculamos las medidas de precisión \(MSE\) y \(RMSE\).
dat.predict |>
summarise(
MSE = mean((atraso_llegada - predictions)^2),
RSME = sqrt(MSE)
) |>
print()## # A tibble: 1 × 2
## MSE RSME
## <dbl> <dbl>
## 1 371. 19.3Tal y como habíamos visto con la técnica de validación, el valor del \(RSME\) obtenido tras el proceso de predicción, es muy alto, de manera que las predicciones obtenidas tienen un error medio de \(\pm\) 19.27 minutos respecto a los valores reales.