library(tidyverse)
library(modelr)
library(nycflights13)
library(lubridate)
library(datos)
library(ggdark)
library(ggrepel)
library(RColorBrewer)
library(patchwork)Tema personalizado.
dark_theme <- list(dark_theme_gray(base_family = "Fira Sans Condensed Light", base_size = 14) +
theme(plot.title = element_text(family = "Fira Sans Condensed", size=rel(0.94), color='grey70'),
plot.background = element_rect(fill = "black"),
plot.caption = element_text(family = "Fira Sans Condensed", size=rel(0.7), color='grey50'),
axis.title.y = element_text(family = "Fira Sans Condensed", size=rel(0.84), color='grey50'),
axis.title.x = element_text(family = "Fira Sans Condensed", size=rel(0.84), color='grey50'),
axis.text.x = element_text(family = "Fira Sans Condensed", size=rel(0.84), color='grey50'),
axis.text.y = element_text(family = "Fira Sans Condensed", size=rel(0.84), color='grey50'),
legend.title = element_text(family = "Fira Sans Condensed", size=rel(0.84), color='grey50'),
legend.text = element_text(family = "Fira Sans Condensed", size=rel(0.8), color='grey50'),
panel.background = element_blank(),
panel.grid.major = element_line(color = "grey30", size = 0.2),
panel.grid.minor = element_line(color = "grey30", size = 0.2),
legend.background = element_blank(),
axis.ticks = element_blank(),
legend.key = element_blank(),
),
scale_color_brewer(palette="Spectral"))
update_geom_defaults("point", list(colour = "deeppink3"))
update_geom_defaults("line", list(colour = "deepskyblue3"))
update_geom_defaults("boxplot", list(colour = "deepskyblue3"))
update_geom_defaults("text_repel", list(box.padding = 0.7,
label.padding = 0,
point.padding = 0.5,
segment.color = 'grey50',
label.size = NA,
fill = "black",
alpha = 0.9))Los datos de vuelos.
Contamos la cantidad de vuelos por día y los graficamos.
vuelos_por_dia <- vuelos %>%
mutate(fecha = make_date(anio, mes, dia)) %>%
group_by(fecha) %>%
summarise(n=n())
head(vuelos_por_dia)## # A tibble: 6 x 2
## fecha n
## <date> <int>
## 1 2013-01-01 842
## 2 2013-01-02 943
## 3 2013-01-03 914
## 4 2013-01-04 915
## 5 2013-01-05 720
## 6 2013-01-06 832ggplot(vuelos_por_dia, aes(x=fecha, y=n)) +
geom_line(size=1) +
dark_theme +
ggtitle('Cantidad de vuelos de NYC por día')
Qué está pasando?
Grafiquemos los vuelos por día de la semana y veamos.
vuelos_por_dia <- vuelos_por_dia %>%
mutate(dia_semana = wday(fecha, label=T))ggplot(vuelos_por_dia, aes(x=dia_semana, y=n)) +
geom_boxplot() +
dark_theme +
ggtitle("Cantidad de vuelos cada dia de la semana.")
Grafiquemos los histogramas para cada dia.
ggplot(vuelos_por_dia, aes(n)) +
geom_histogram(color='white', fill='deeppink4') +
facet_wrap(~ dia_semana) +
dark_theme## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Superpongamos los histogramas.
ggplot(vuelos_por_dia, aes(n)) +
geom_histogram(aes(color=dia_semana, fill=dia_semana), alpha = 0.3) +
dark_theme +
theme(legend.position = 'right')## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Hay un fuerte efecto del día de la semana. En particular los fines de
semana cae el número de vuelos. ¿Por qué? Una hipótesis: La mayoría son viajes de negocios.
Tratamos de remover el patrón, utilizando un modelo.
mod <- lm(n ~ dia_semana, data=vuelos_por_dia)
cuadricula <- vuelos_por_dia %>%
data_grid(dia_semana) %>%
add_predictions(mod, "n")Graficamos el modelo:
ggplot(vuelos_por_dia, aes(dia_semana, n)) +
geom_point(color='deepskyblue4') +
geom_boxplot(alpha=0.1) +
geom_point(data = cuadricula, size= 4) +
dark_theme +
ggtitle('Predicción del modelo a los dias de la semana.')
# Obs. El modelo es suceptible a valores atípicos.
Que tan bueno resultó el modelo?
Veamos los residuos:
Agregamos los residuos al dataframe.
vuelos_por_dia <- vuelos_por_dia %>%
add_residuals(mod)Los graficamos:
ggplot(vuelos_por_dia, aes(x=fecha, y=resid)) +
geom_ref_line(h=0, colour='grey50', size=0.5) +
geom_point() +
geom_line(alpha=0.6) +
dark_theme +
ggtitle("Residuos del modelo lineal")
# Que está pasando? Parece que hay algunas fechas donde la cantidad de vuelos disminuye drásticamente. # Separemos por dias.
ggplot(vuelos_por_dia, aes(fecha, resid, color=dia_semana)) +
geom_ref_line(h=0, colour='grey50', size=0.5) +
geom_line(alpha=0.6) +
geom_point() +
dark_theme 
Obs:
1) Los sábados parecen seguir otra tendencia. El modelo sobreestima significativamente (>100) la cantidad de vuelos en algunos dias,
en particular los viernes y domingos.
2) El modelo comete errores sistemáticos en el verano (Junio, Julio y Agosto), donde subestima la cantidad de vuelos en estas fechas.
Se puede mejorar?
Filtermos los dias anómalos donde el modelo sobreestima los vuelos y veamos si podemos encontrar algún patrón.
dias_raros <- vuelos_por_dia %>% filter(abs(resid) > 100)Agregemos a la gráfica.
p <- ggplot(vuelos_por_dia, aes(fecha, resid, color=dia_semana)) +
geom_ref_line(h=0, colour='grey50', size=0.5) +
geom_line(alpha=0.6) +
geom_point() +
dark_theme
p + geom_point(data=dias_raros, aes(fecha, resid), size=4, show.legend=F)+
geom_label_repel(data=dias_raros, aes(label=format(fecha, format='%e %b')), show.legend=F) +
dark_theme +
ggtitle("Fechas anómalas")## Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
# Ejercicio 1: ¿Por qué hay menos vuelos el 20 de enero, el 26 de mayo y el 1ero de septiembre? Se generaliza para otros años? # Son feriados no recurrentes, por ejemplo el 20 de enero se ignaguró el segudno mandato de Obama.
Las otras fechas parecen ser feriados recurrentes.
Hagamos un suavizado para reslatar la tendencia que hablamos.
vuelos_por_dia %>%
ggplot(aes(fecha, resid)) +
geom_ref_line(h=0, colour='grey50', size=0.5) +
geom_line() +
geom_smooth(se=F, span=0.3, method = "loess") +
dark_theme +
ggtitle("Tendencia anómala")## `geom_smooth()` using formula 'y ~ x'
Veamos que occure los sábados y el posible efecto estacional sobre la cantidad de vuelos.
vuelos_por_dia %>%
filter(dia_semana == 'Sat') %>%
ggplot(aes(fecha, resid)) +
geom_ref_line(h=0, colour='grey50', size=0.5) +
geom_line(alpha=1, size=1) +
geom_point()+
dark_theme +
scale_x_date(date_breaks = "1 month", date_labels = '%b') +
ggtitle("Cantidad de vuelos los Sábados")
## Obs. # Una posible explicación de porqué la gente tiende a viajar más los sábados en el verano, # es que se van de vacaciones. Las vacaciones de verano en el 2013 en NYC fueron del 26 de junio hasta el # 9 de setiembre.
Tambien podemos ver que hay más vuelos en el perido de primavera (Febrero-Mayo) que en otoño (Setiembre-Diciembre).
Probemos esta hipótesis. Construimos una variable categórica para las estaciones y vemos si nos ayuda a modelar los datos.
trimestre <- function(fecha){
cut(fecha,
breaks = ymd(20130101, 20130605, 20130825, 20140101),
labels = c("primavera", "verano", "otoño")
)
}
vuelos_por_dia <- vuelos_por_dia %>%
mutate(trimestre = trimestre(fecha))Filtremos los sábados por temporada y grafiquemos.
vuelos_por_dia %>%
filter(dia_semana == "Sat") %>%
ggplot(aes(fecha, resid, color = trimestre)) +
geom_point() +
geom_line()+
scale_x_date(date_breaks="1 month", date_labels = "%b") +
dark_theme +
ggtitle("Vuelos los sábados por temporada.")
# Que ocurre con los otros días de semana?
vuelos_por_dia %>%
filter(dia_semana == "Fri") %>%
ggplot(aes(fecha, resid, color = trimestre)) +
geom_point() +
geom_line()+
scale_x_date(date_breaks="1 month", date_labels = "%b") +
dark_theme +
ggtitle("Vuelos los viernes por temporada.")
# Los fines de semana:
vuelos_por_dia %>%
filter(dia_semana %in% c("Fri","Sun","Sat")) %>%
ggplot(aes(n)) +
geom_histogram(aes(color=dia_semana, fill=dia_semana), alpha = 0.3) +
dark_theme +
theme(legend.position = 'right') +
facet_wrap(~trimestre)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Todos los dias:
vuelos_por_dia %>%
ggplot(aes(dia_semana, n, color = trimestre)) +
geom_boxplot() +
geom_point(data=cuadricula, color = 'deeppink3', shape=95, size=12) +
dark_theme
El modelo previo parece describir mejor los datos de la primavera.
Probemos con otros modelo, utilizando la variable trimestre.
mod1 <- lm(n ~ dia_semana, data = vuelos_por_dia)
mod2 <- lm(n ~ dia_semana*trimestre, data = vuelos_por_dia)Comparemos los residuos de los modelos.
vuelos_por_dia %>%
gather_residuals(sin_trimestre = mod1, con_trimestre=mod2) %>%
ggplot(aes(fecha, resid, color = model)) +
geom_line(alpha=0.7, size=1) +
geom_point() +
dark_theme +
scale_color_manual(values=c("cyan3", "deeppink3")) ## Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
No parece que ajusta mejor.
Comparemos los modelos en los datos.
cuadricula <- vuelos_por_dia %>%
data_grid(dia_semana, trimestre) %>%
spread_predictions(mod1, mod2)
ggplot(vuelos_por_dia, aes(dia_semana, n, color=trimestre)) +
geom_boxplot()+
geom_point(data=cuadricula, aes(x=dia_semana, y=mod1), color="deeppink3", shape=95, size=10)+
geom_point(data=cuadricula, aes(x= dia_semana, y=mod2), color = "cyan3", shape=95, size=10)+
dark_theme +
theme(legend.position = "none") +
facet_wrap(~ trimestre)
El modelo nuevo ajusta mejor al valor de las medias de los datos, en particular en primavera y verano.
Sin embargo la preponderancia de datos atípicos en otoño lleva a que ambos modelos tengan residuos significativos.
Probemos crear un modelo más robusto a la presencia de datos atípicos.
mod3 <- MASS::rlm(n ~ dia_semana*trimestre, data = vuelos_por_dia)
vuelos_por_dia %>%
gather_residuals(sin_trimestre = mod1, con_trimestre=mod2, robusto=mod3) %>%
ggplot(aes(fecha, resid, color = model)) +
geom_line(alpha=0.7, size=1) +
geom_point() +
dark_theme +
scale_color_manual(values=c("deeppink3", "green3", "cyan3")) ## Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Bastante mejor, en particular en el verano donde se logró reducir el error sistemático observado en ésta temporada.
Otro enfoque: Podemos aumentar la capacidad del modelo, abandonar los modelos lineales y por ejemplo utilizar splines.
Generalmente ganamos predictibilidad a costo de interpretabilidad.
library(splines)
spline_mod <- MASS::rlm(n ~ dia_semana*ns(fecha, 5), data = vuelos_por_dia)
vuelos_por_dia %>%
data_grid(dia_semana, fecha = seq_range(fecha, n = 13)) %>%
add_predictions(spline_mod) %>%
ggplot(aes(fecha, pred, colour = dia_semana)) +
geom_line() +
geom_point() +
dark_theme
# Nuevamente aparece el patrón de los sábados. # Comparemos el modelo de splines con los otros.
vuelos_por_dia %>%
gather_residuals(sin_trimestre = mod1, con_trimestre=mod2, robusto=mod3, splines = spline_mod) %>%
ggplot(aes(fecha, resid, color = model)) +
geom_line(alpha=0.7, size=1) +
geom_point() +
dark_theme +
scale_color_manual(values=c("deeppink3", "green3", "cyan3", "gold2")) ## Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
Logra captar el efecto de temporada en el verano, pero aparecen nuevos problemas en invierno.
cuadricula <- vuelos_por_dia %>%
data_grid(dia_semana, fecha=seq_range(fecha, n = 1), trimestre) %>%
spread_predictions(mod1, mod2, mod3, spline_mod)
ggplot(vuelos_por_dia, aes(dia_semana, n, color=trimestre)) +
geom_boxplot()+
geom_point(data=cuadricula, aes(x=dia_semana, y=mod1), color="deeppink3", shape=95, size=10)+
geom_point(data=cuadricula, aes(x= dia_semana, y=mod2), color = "cyan3", shape=95, size=10)+
geom_point(data=cuadricula, aes(x=dia_semana, y=mod3), color="green3", shape=95, size=10)+
geom_point(data=cuadricula, aes(x=dia_semana, y=spline_mod), color="gold2", shape=95, size=10)+
dark_theme +
theme(legend.position = "none") +
facet_wrap(~ trimestre)
No aproxima bien a la media de los datos, parece que los datos atípicos están causando errores, habría que refinar el modelo.
Ejercicio 2:
Que representan esos 3 días con altos residuos positivos? Se generaliza a otros años?
vuelos_por_dia %>%
slice_max(n=3, resid) %>%
mutate(descr_fecha = format(fecha, format = "%e %b"))## # A tibble: 3 x 6
## fecha n dia_semana resid trimestre descr_fecha
## <date> <int> <ord> <dbl> <fct> <chr>
## 1 2013-11-30 857 Sat 112. otoño "30 Nov"
## 2 2013-12-01 987 Sun 95.5 otoño " 1 Dec"
## 3 2013-12-28 814 Sat 69.4 otoño "28 Dec"El 30 de Nov, 1 Dec y 28 de Dic son días despúes de feriados recurrentes, se generaliza a otros años.
Ejercico 3:
Crea una nueva variable que divida dia_semana en periodos pero sólo para sábados. O sea tener: Thu, Fri, Sat-verano, Sat-primavera, Sat-otoño.
Cómo se compara con los otros?
trimestre <- function(fecha){
cut(fecha,
breaks = ymd(20130101, 20130605, 20130825, 20140101),
labels = c("primavera", "verano", "otoño")
)
}
dias_tri <- vuelos_por_dia %>%
mutate(trimestre = trimestre(fecha)) %>%
mutate(dias = ifelse(dia_semana == "Sat", paste0(dia_semana, "-", trimestre), as.character(dia_semana)))
head(dias_tri, 14)## # A tibble: 14 x 6
## fecha n dia_semana resid trimestre dias
## <date> <int> <ord> <dbl> <fct> <chr>
## 1 2013-01-01 842 Tue -109. primavera Tue
## 2 2013-01-02 943 Wed -19.7 primavera Wed
## 3 2013-01-03 914 Thu -51.8 primavera Thu
## 4 2013-01-04 915 Fri -52.5 primavera Fri
## 5 2013-01-05 720 Sat -24.6 primavera Sat-primavera
## 6 2013-01-06 832 Sun -59.5 primavera Sun
## 7 2013-01-07 933 Mon -41.8 primavera Mon
## 8 2013-01-08 899 Tue -52.4 primavera Tue
## 9 2013-01-09 902 Wed -60.7 primavera Wed
## 10 2013-01-10 932 Thu -33.8 primavera Thu
## 11 2013-01-11 930 Fri -37.5 primavera Fri
## 12 2013-01-12 690 Sat -54.6 primavera Sat-primavera
## 13 2013-01-13 828 Sun -63.5 primavera Sun
## 14 2013-01-14 928 Mon -46.8 primavera MonEl modelo:
mod_sab <- lm(n ~ dias, data = dias_tri)
mod1 <- lm(n ~ dia_semana, data = dias_tri)
mod2 <- lm(n ~ dia_semana*trimestre, data = dias_tri)
dias_tri %>%
gather_residuals(sin_trimestre = mod1, con_trimestre=mod2, sabados=mod_sab) %>%
ggplot(aes(fecha, resid, color = model)) +
geom_line(alpha=0.7, size=1) +
geom_point() +
dark_theme +
scale_color_manual(values=c("deeppink3", "green3", "cyan3")) ## Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
# No hay mejora substancial.
Ejercicio 4:
Crea una variable dia_semana que combina el dia de la semana, periodos (para sábados), y feriados públicos. ¿Cómo se ven los residuos?
dias_feriados <- dias_tri %>%
mutate(feriados = case_when(fecha %in% ymd(c(20130101, 20130121, 20130218, 20130527, 20130704,
20130902, 20131028, 20131111, 20131128, 20131225)) ~ "feriado", TRUE ~ "None")) %>%
unite(trimestre, dias, feriados)El modelo:
mod_feriados <- lm(n ~ dia_semana*trimestre, data = dias_feriados)
mod1 <- lm(n ~ dia_semana, data = dias_feriados)
dias_feriados %>%
gather_residuals(feriados=mod_feriados, base=mod1) %>%
ggplot(aes(fecha, resid, color=model)) +
geom_line() +
geom_point() +
dark_theme
# Hay una leve mejora, casi imperceptible.
Ejercicio 6:
Que sucede si ajustas un efecto de día de la semana que varía según el mes o varía mes a mes (es decir n ~ dia_semana*month)? ¿Por qué no es útil?
vuelos_mes <- vuelos_por_dia %>%
mutate(mes = month(fecha))
mod1 <- lm(n ~ dia_semana*trimestre, data=vuelos_mes)
mod_mes <- lm(n ~ dia_semana*mes, data= vuelos_mes)
vuelos_mes %>%
gather_residuals(base = mod1,mes = mod_mes) %>%
ggplot(aes(fecha, resid, color=model)) +
geom_line() +
geom_point() +
dark_theme
# Es peor. Considera ventanas más grandes (meses) y termina con más cantidad de valores atípicos.
Ejercicio 7.
Una hipótesis era que las personas que salen los domingos, viajan por negocio, ¿hay evidencia de esto?
Si es verdad se esperaría ver más vuelos en la tarde del domigno a lugares lejanos.
domingos <- vuelos %>%
filter(hora > 19, distancia > mean(distancia))
vuelos_domingos <- domingos %>%
mutate(fecha = make_date(anio, mes, dia)) %>%
group_by(fecha) %>%
summarise(n=n())
vuelos_domingos <- vuelos_domingos %>%
mutate(dia_semana = wday(fecha, label=T))
p1 <- ggplot(vuelos_por_dia, aes(x=dia_semana, y=n)) +
geom_boxplot() +
dark_theme
p2 <- ggplot(data=vuelos_domingos, aes(x=dia_semana, y = n)) +
geom_boxplot() +
dark_theme
p1 + p2