```{r}
library(tidyverse)
library(modelr)
library(nycflights13)
library(lubridate)
library(datos)
library(ggdark)
library(ggrepel)
library(RColorBrewer)
library(patchwork)
```

## Tema personalizado.
```{r}
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.

```{r}
vuelos_por_dia  <- vuelos  %>% 
    mutate(fecha = make_date(anio, mes, dia))  %>% 
    group_by(fecha)  %>% 
    summarise(n=n())

head(vuelos_por_dia)
```

```{r}
ggplot(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.

```{r}
vuelos_por_dia  <- vuelos_por_dia  %>% 
    mutate(dia_semana = wday(fecha, label=T))
```

```{r}
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.
```{r}
ggplot(vuelos_por_dia, aes(n)) +
    geom_histogram(color='white', fill='deeppink4') +
    facet_wrap(~ dia_semana) +
    dark_theme
```

# Superpongamos los histogramas.
```{r}
ggplot(vuelos_por_dia, aes(n)) +
    geom_histogram(aes(color=dia_semana, fill=dia_semana), alpha = 0.3) +
    dark_theme +
    theme(legend.position = 'right')
```

# 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.
```{r}
mod  <-  lm(n ~ dia_semana, data=vuelos_por_dia)    

cuadricula  <-  vuelos_por_dia  %>% 
    data_grid(dia_semana)  %>% 
    add_predictions(mod, "n")
```

# Graficamos el modelo:
```{r}
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.
```{r}
vuelos_por_dia  <- vuelos_por_dia  %>% 
    add_residuals(mod)
```

# Los graficamos:
```{r}
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.
```{r}
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.
```{r}
dias_raros <- vuelos_por_dia  %>% filter(abs(resid) > 100)
```

# Agregemos a la gráfica.
```{r}
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")
```
# 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.
```{r}
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")
```

# Veamos que occure los sábados y el posible efecto estacional sobre la cantidad de vuelos.
```{r}
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.
```{r}
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.
```{r}
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?
```{r}
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:
```{r}
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)
```
# Todos los dias:
```{r}
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.

```{r}
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.
```{r}
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")) 
```

# No parece que ajusta mejor. 
# Comparemos los modelos en los datos.

```{r}
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.
```{r}
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")) 
```

# 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.

```{r}
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.

```{r}
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")) 
```

# Logra captar el efecto de temporada en el verano, pero aparecen nuevos problemas en invierno.

```{r}
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?

```{r}
vuelos_por_dia  %>% 
    slice_max(n=3, resid)  %>% 
    mutate(descr_fecha = format(fecha, format = "%e %b"))
```
# 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?

```{r}
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)
```
# El modelo:
```{r}
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")) 
```
# 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?
```{r}
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:
```{r}
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?
```{r}
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.
```{r}
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
```