El problema de los vuelos diarios

pacman::p_load(tidyverse, modelr, nycflights13, lubridate)

Business question: ¿qué afecta al número de vuelos de un día?

Empecemos para hacer un filtrado para saber el número de vuelos que salen en un día:

daily <- flights %>%
  mutate(date = make_date(year, month, day)) %>% # 
  group_by(date) %>%
  summarise(n = n()) # creamos el contador
head(daily, 10)

Visualizamos la información:

daily %>%
  ggplot(aes(date, n)) + geom_line()

Creación de la hipótesis

Sacamos una hipótesis; puede que haya una relación del número de vuelos con el día de la semana, viendo que el lunes sube y el domingo bajan. Por lo tanto puede que el día de la semana sea la variable más importante.

Creamos una columna adicional con el día de la semana:

daily %>%
  mutate(wday = wday(date, label = TRUE, locale = "us")) -> daily # creamos el numero de la semana
head(daily)

Visualizamos el día de la semana con el conteo:

daily %>%
  ggplot(aes(wday, n)) +
    geom_boxplot()

Conclusiones de la visualización:

Hay menos vuelos los sabados y los domingos, dándonos a entender que la mayoría de persona que pasan por el aeropuerto de Nueva York es por trabajo.

Creación del modelo

Vamos a intentar predecir el número de vuelos que salen solo en base al día de la semana:

# Creamos el modelo
mod1 <- lm(n ~ wday, data = daily)
# Creación del grid para comprovar resultados
grid <- daily %>%
  data_grid(wday) %>%
  add_predictions(mod1, "n")
# Visualización de los resultados
daily %>%
  ggplot(aes(wday, n)) +
    geom_boxplot() +
    geom_point(data = grid, color = "red", size = 4)

La predicción esta marcado en rojo, y es lo que nos valora el modelo.

Análisis de residuos

# Añadimos los residuos
daily <- daily %>%
  add_residuals(mod1)
# Visualizamos los residuos
daily %>%
  ggplot(aes(date, resid)) +
    geom_ref_line(h = 0) + # añadimod la referencia en 0 para valorar la información
    geom_line()

Viendo los errores del modelo en relación a los días de la semana, vemos que los picos se han reducido respecto al gráfico de lineas anterior, pero siguen habiendo picos. Puede que el modelo siga sin explicar determinados días de la semana. Hay patrones que siguen apareciendo.

En junio aparece el primer pico
En septiembre y octubre se acentuan
En navidades hay el mismo problema

Dividamos el gráfico en día de la semana por colores para obtener una mayor información:

daily %>%
  ggplot(aes(date, resid, color = wday)) +
    geom_ref_line(h = 0) + 
    geom_line()

Podemos ver que nuestro modelo falla:

En los sábados, concretamente en verano, donde están muy infravalorados, y en septiembre y octubre están sobrevalorados. Se debe mejorar el factor estacionalidad
Hay días que tienen muchos menos vuelos que esperávamos:

daily %>%
  filter(abs(resid) > 100)

Aquí podemos ver rápidamente porqué algunos días tienen un error muy grande:

2013-07-04: el 4 de junio es el día más importante de estados unidos, es normal que no lo pueda predecir.
Acción de gracias también aparece.
Navidad, nochebuena y fin de año […]

Podemos ver que los picos son factores sociopolíticos, impredecibles para el modelo. Pero podríamos añadir alguna restricción adicional. Pero hay una solución, crear una versión suavizado de la tendencia para todo el año con el objetivo de reducir el factor social:

daily %>%
  ggplot(aes(date, resid)) +
    geom_ref_line(h = 0) +
    geom_line(color = "grey30") +
    geom_smooth(se = TRUE, span = 0.2)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Conclusiones

En invierno hay menos vuelos de los que el modelo espera

Vamos a mejorar el modelo porque nos estamos dejando información.

Vuelos por día de la semana y trimestre

Vamos a intentar mejorar la eficacia únicamente los sábados:

daily %>%
  filter(wday == "Sat") %>%
  ggplot(aes(date, n)) +
    geom_point() +
    geom_line() +
    scale_x_date(NULL, date_breaks = "1 month", date_labels = "%b") # con este parámetro añadiremos un cambio en el factor de escala x en formato  fecha,donde no añadimos ningún nombre adicional, el date breaks quiero que sea en grupos de un mes, y el date breaks quiero que sea el día del mes versión reducida por lo tanto uso la sintaxis de lubridate

Usamos el formato de puntos y lineas para especificar es una interpoblación (sólo sábados).

Sospechamos que el factor que provoca que nuestro modelo falle son las vacaciones de verano, desde junio hasta septiembre (26 de junio hasta el 9 de septiembre).

También parece que haya más vuelos un sábado en primavera que un sabado en otoño. Hipótesis; puede que la gente se reserve en otoño para volar en navidades.

Análisis vacaciones de 5 de Junio al 9 de septiembre

Usaremos desde el 5 de junio, la primera semana de junio, hasta el 25 de agosto, que es cuando empieza a bajar el número de vuelos.

term <- function(date) { # dada una fecha cualquiera
  cut(date,
      breaks = ymd(20130101,20130605,20130825,20140101),
      labels = c("spring", "summer", "fall")
  )
}

daily <- daily %>%
  mutate(term = term(date))

daily %>%
  filter(wday == "Sat") %>%
  ggplot(aes(date, n, color = term))+
    geom_point(alpha = .25) +
    geom_line() +
    scale_x_date(NULL, date_breaks = "1 month", date_labels = "%b")

En esta gráfica vemos claramente como el verano tiene muchos vuelo, en otoño muy bajos, a pesar de que hayan las vacaciones de invierno, y primavera estaría en el medio.

Vamos a ver como estas varaibles afectan a los trimestres:

daily %>%
  ggplot(aes(wday, n, color = term)) +
  geom_boxplot()

Las cajas de verano son más altas que el resto. Esta variancia significativa a lo largo de los trimestras nos dice que la predicción del número de vuelos no puede ir solo son el día de la semana, sinó que debe ir también acompañada del trimestre en que estemos.

Creación del segundo modelo

Vamos a ver como interactua los trimestres con el número de vuelos:

mod2 <- lm(n ~ wday*term, data = daily)

Interpretación de resudios

daily  %>%
  gather_residuals(without_term = mod1, with_term = mod2) %>%
  ggplot(aes(date, resid, color = model)) +
    geom_line(alpha = 0.7)

El modelo con trimestre esta más cerca del cero en verano, así que parece que ha mejorado considerablemente. Seguimos con el análisis:

# Creamos la parrilla para hacer las predicciones
grid <- daily %>%
  data_grid(wday, term) %>%
  add_predictions(mod2, "n") 

daily %>%
  ggplot(aes(wday, n)) +
    geom_boxplot() +
    geom_point(data = grid, color = "red") +
    facet_wrap(~term)

En primavera no tenemos outliers, pero aún siguen habiendo muchos outliers que provocan a nuestras previsiones desplazarse mucho del centro de las cajas. Esto es debido a que hemos usado un modelo lineal. Con el objetivo de reducir su efecto, usaremos otro modelo más robusto.

Uso de modelo lineal robusto

Usaremos este sistema para reducir el efecto de los outliers. Usaremos el paquete MASS

mod3 <- MASS::rlm(n ~ wday * term, data = daily)

daily %>%
  add_residuals(mod3, "resid") %>%
  ggplot(aes(date, resid)) +
    geom_hline(yintercept = 0, size = 2, color = "white") +
    geom_line()

Podemos ver como hemos aplanado mucho la tendencia, en verano es aproximadamente 0. Aun así siguen habiendo muchos outlier en invierno.

Modelos múltiples

Es muy buena práctica computar variables siempre que vamos trabajando, es decir, ir creando un conjunto de variables dentro de una función para no equivocarse.

# Creamos una función para ir guardando las variables que creamos que son representativas del modelo
compute_vars <- function(data) {
  data %>%
    mutate(
      term = term(date),
      wday = wday(date, labels = TRUE)
    )
}

# otra opción es poner la transformación directamente (sin simplificar) en el propio modelo
wday2 <- function(x) wday(x, labels = TRUE)
# mod4 <- lm(n ~ wday2(date) * term(date), data = daily)

Vamos a ver si podemos usar el día de la semana con la fecha con el objetivo de tratar el año al completo:

library(splines)
mod <- MASS::rlm(n ~ wday * ns(date, 5), data = daily) # usaremos el orden 5

daily %>%
  data_grid(wday, date = seq_range(date, n = 13)) %>% # cogeremos 13 valores, uno para cada día del año, cogiendo así como último el uno de enero de 2014. Será una  parrilla de datos con 13  fechas y 7 días de a semana
  add_predictions(mod) %>%
  ggplot(aes(date, pred, color = wday)) +
    geom_line() +
    geom_point()

Casos Reales_Modelos con el dataset Vuelos

Joan Claverol _ Curso Juan Gabriel Gomila _ Udemy

15 de agosto de 2018