1、计算每日航班数量

计算每日航班数量

> daily <-  flights %>%
+   mutate(date=make_date(year,month,day)) %>%
+   group_by(date) %>%
+   summarize(n=n())
`summarise()` ungrouping output (override with `.groups` argument)
> daily
# 可视化
> ggplot(daily,aes(date,n))+
+   geom_line()

image.png

• n()的用法:分组计数

> daily <-  flights %>%
+   mutate(date=make_date(year,month,day)) %>%
+   group_by(date)
> daily
> daily
# A tibble: 336,776 x 20
# Groups:   date [365]
    year month   day dep_time sched_dep_time dep_delay arr_time sched_arr_time
   <int> <int> <int>    <int>          <int>     <dbl>    <int>          <int>
 1  2013     1     1      517            515         2      830            819
 2  2013     1     1      533            529         4      850            830
 3  2013     1     1      542            540         2      923            850
 4  2013     1     1      544            545        -1     1004           1022
 5  2013     1     1      554            600        -6      812            837
 6  2013     1     1      554            558        -4      740            728
 7  2013     1     1      555            600        -5      913            854
 8  2013     1     1      557            600        -3      709            723
 9  2013     1     1      557            600        -3      838            846
10  2013     1     1      558            600        -2      753            745
> count(group_by(daily,date))
# A tibble: 365 x 2
# Groups:   date [365]
   date           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   832
 7 2013-01-07   933
 8 2013-01-08   899
 9 2013-01-09   902
10 2013-01-10   932

数据显然有以周为单位的变化，这影响了数据的长期观察。

检查航班数量在每一天正宗的分布

> Sys.setlocale("LC_ALL","English")
[1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"
> daily <-  daily %>%
+   mutate(wday=wday(date,label = T,locale = Sys.getlocale(category = "LC_TIME")))
> daily
# A tibble: 365 x 3
   date           n wday 
   <date>     <int> <ord>
 1 2013-01-01   842 Tue  
 2 2013-01-02   943 Wed  
 3 2013-01-03   914 Thu  
 4 2013-01-04   915 Fri  
 5 2013-01-05   720 Sat  
 6 2013-01-06   832 Sun  
 7 2013-01-07   933 Mon  
 8 2013-01-08   899 Tue  
 9 2013-01-09   902 Wed  
10 2013-01-10   932 Thu  
# ... with 355 more rows
> ggplot(daily,aes(wday,n))+
+   geom_boxplot()

image.png

必须把locale设为英文环境，否则lucbridate会根据当前环境读取中文，显示中文的星期几

消除强烈的模式——建立模型

> ## 拟合模型
> mod <-  lm(n~wday,data = daily)
> grid <-  daily %>% 
+   data_grid(wday) %>% 
+   add_predictions(mod,"n")
> ggplot(daily,aes(wday,n))+
+   geom_boxplot()+
+   geom_point(data = grid,aes(wday,n),color="red",size=4)
> ggplot(daily,aes(wday,n))+
+   geom_boxplot()+
+   geom_point(data = grid,color="red",size=4)
> grid <-  daily %>% 
+   data_grid(wday) %>% 
+   add_predictions(mod,"n")
> ggplot(daily,aes(wday,n))+
+   geom_boxplot()+
+   geom_point(data = grid,color="red",size=4)

image.png

分析残差

> daily <- daily %>% 
+   add_residuals(mod) 
> daily %>% 
+   ggplot(aes(date,resid))+
+   geom_ref_line(h=0)+
+   geom_line()

image.png

六月开始模型并不适用
→进一步分析——按照wday分别展示

按照wday分别展示

daily %>% 
  ggplot(aes(date,resid,color=wday))+
  geom_ref_line(h=0)+
  geom_line()
# 显示航班特别少的日期
daily %>% 
  filter(resid < -100)

image.png

显示更平滑的长期趋势: smooth()

daily %>% 
  ggplot(aes(date,resid))+
  geom_ref_line(h=0)+
  geom_line(color="grey50")+
  geom_smooth(se=F,span=0.20)

image.png

季节性星期六效应

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

image.png

labels = b% : 见strptime {base}
%b
Abbreviated month name in the current locale on this platform. (Also matches full name on input: in some locales there are no abbreviations of names.)

从图中可看出周六航班的季节变化。可能与学期假期相关。

按学期分类处理数据

> term <-  function(date){
+   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=1/3)+
+   geom_line()+
+   scale_x_date(
+     NULL,
+     date_breaks = "1 month",
+     date_labels = "%b")

image.png

查看学期变量如何影响一周中其他wday

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

image.png

拟合去除每学期周内效应的模型

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

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

image.png

将预测值覆盖到数据上

grid <-  daily %>% 
  data_grid(wday,term) %>% 
  add_predictions(mod2,"n")

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

image.png

发现问题：离群点
解决：使用MASS::rlm()

处理离群点—— MASS::rlm()

library(MASS)
mod3 <- 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()

image.png