R数据科学-第十二章模型构建
在本章中, 我们将基于真实的数据显示, 详细阐述从模型 bin 到数据的理解过程
首先运用数据可视化方法识别出模型。接着通过该模型更细致且准确地提取出特征模式。最后进一步分析剩余的误差项。
为了便于分析ggplot2::diamonds和nycflights13::flights数据集的特征关系, 所需的R包如下:
library(tidyverse)
library(modelr)
library(ggplot2)
options(na.action=na.warn)
library(nycflights13)
library(lubridate)一、为什么质量差的钻石更贵
通过数据可视化技术我们揭示出：外观不佳色调不匀或杂质含量高的钻石其定价水平显著高于其他 diamonds with better cut polish and higher clarity.
> ggplot(diamonds,aes(cut,price))+geom_boxplot()
> ggplot(diamonds,aes(color,price))+geom_boxplot()D到J颜色等级降低。
> ggplot(diamonds,aes(clarity,price))+geom_boxplot()I1纯度最差。
造成这个反常现象的原因是由于一个混淆变量的存在：重量。
重量是决定钻石价格最重要的因素，而质量差的钻石往往更重些。
> ggplot(diamonds,aes(carat,price))+geom_point()
利用一个模型来提取carat变量的作用, 从而使得其他变量的影响更加明显
第一步：数据预处理
1.重点关注carat<2.5的数据，因为那占了绝大部分数据。
2.将carat和price数据进行对数转换，以便更容易看出二者之间的关系。
> diamonds2<-diamonds%>%filter(carat<=2.5)%>%mutate(lprice=log2(price),lcarat=log2(carat))
> ggplot(diamonds2,aes(lcarat,lprice))+geom_point()
第二步：构建模型去除二者之间这种强烈的线性模式
1.拟合一个模型让这种模式变为显式的
#首先基于原始carat的最大最小值生成均匀分布的20个数据点，并在这些数据中添加lcarat=log2(carat)这一列特征。接着利用该新增特征与预训练模型mod_diamonds进行拟合，并输出预测结果命名为lprice。最后通过对lprice进行指数变换完成最终预测值的还原。
> grid<-diamonds2%>%data_grid(carat=seq_range(carat,20))%>%mutate(lcarat=log2(carat))%>%add_predictions(mod_diamonds,"lprice")%>%mutate(price=2^lprice)
> ggplot(diamonds2,aes(carat,price))+geom_point()+geom_line(data=grid,color='red',size=1)
2.检查残差，验证模型的好坏：
> diamonds2<-diamonds2%>%add_residuals(mod_diamonds,"lresid")
> ggplot(diamonds2,aes(lcarat,lresid))+geom_point()
3.用残差替换price，以移除线性关系：
由于价格可被分解为与carat相关的部分及残差部分两个方面，在这种情况下采用残差替代原始价格，则等同于直接排除了与carat相关的那一部分
> ggplot(diamonds2,aes(cut,lresid))+geom_boxplot()
> ggplot(diamonds2,aes(color,lresid))+geom_boxplot()
> ggplot(diamonds2,aes(clarity,lresid))+geom_boxplot()
注意：我们估计到的模型实际上为：
当我们用残差去代替因变量时实际上代替的是log2(price)。
故图形中，残差=-1代表log2(price)=-1,此时price=1/2，即该点价格为预计价格的一半。
二、哪些因素影响了每日航班数量
第一步：计算每天出发的航班数量
> daily<-flights%>%mutate(date=make_date(year,month,day))%>%group_by(date)%>%summarize(n=n())
> daily
# A tibble: 365 x 2
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
# ... with 355 more rows可视化
> ggplot(daily,aes(date,n))+geom_line()
第二步：航班数量在每一周的分布
掌握长期趋势并非易事,因为数据呈现明显的日间波动特征;因此首先统计各周工作日与休息日的航班数量变化情况；
> daily<-daily%>%mutate(wday=wday(date,label=TRUE))
> ggplot(daily,aes(wday,n))+geom_boxplot()
从图形中可以看出：周末航班的数量相对较少，在这些行程中大部分是用于公差的。这种现象在周六尤其明显，在这种情况下人们通常希望在周一前完成工作安排
第三步：建立模型
> 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,color='red',size=4)
第四步：计算残差并可视化
> daily<-daily%>%add_residuals(mod)
> daily%>%ggplot(aes(date,resid))+geom_ref_line(h=0)+geom_line()
日期与残差的图像已去除绝大多数周内效应。该图片的y轴表示预测值与真实值之间的差异程度。
下面按照wday绘制不同的折线图：
> ggplot(daily,aes(date,resid,color=wday))+geom_ref_line(h=0)+geom_line()
夏季的航班数量要比我们预计的多，秋季的航班数量要比我们预计的少。
下面看一些预计偏差很大的值：
> daily%>%filter(resid<(-100))
# A tibble: 11 x 4
date n wday resid
<date> <int> <ord> <dbl>
1 2013-01-01 842 周二 -109.
2 2013-01-20 786 周日 -105.
3 2013-05-26 729 周日 -162.
4 2013-07-04 737 周四 -229.
5 2013-07-05 822 周五 -145.
6 2013-09-01 718 周日 -173.
7 2013-11-28 634 周四 -332.
8 2013-11-29 661 周五 -306.
9 2013-12-24 761 周二 -190.
10 2013-12-25 719 周三 -244.
11 2013-12-31 776 周二 -175.在过去几天里，航班数量明显低于预期。可以看出这些日期涵盖了美国新年假期（1月1日至1月3日）、独立日（7月4日）、感恩节（11月22日至24日）以及圣诞季（12月初）。
在观察全年数据时似乎呈现出一种更为平滑的趋势，并且通过 geom_smooth 函数能够突出显示这种趋势。
> daily%>%ggplot(aes(date,resid))+geom_ref_line(h=0)+geom_line(color="grey50")+geom_smooth()
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
可以大致看出，5-9月航班较多，1月和12月航班数量较少。
五、季节性星期六效应
我们只关注星期六的航班状况：
> daily%>%filter(wday=="周六")%>%ggplot(aes(date,n))+geom_point()+geom_line()
我们可以观察到六七月间周末航班数量显著增加；这主要归因于暑假时间安排；而不管今天的星期几
春季周六的飞行次数要多于秋季周六的飞行次数可能是由于美国人偏好选择不同时间段出行的原因其中一项原因是人们宁愿选择在感恩节和圣诞节长假期间出游因为这个时间段虽然同样是周末但与秋季其他假期相比更为集中
下面创建一个“学期”变量来粗略地表示学校的三个学期：
> 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=="周六")%>%ggplot(aes(date,n,color=term))+geom_point()+geom_line()
看一下学期因素是如何影响一周中其他天的：
> daily%>%ggplot(aes(wday,n,color=term))+geom_boxplot()
经观察分析可知学期因素对航班数量的影响非常显著；基于此考虑建立一个模型来同时提取学期效应和工作日效应。
> 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()
可以看到效果并没有好太多。
将预测值覆盖在原数据上：
> 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)
由于该模型旨在计算"平均效应"这一指标，在面对存在大量离群点的数据时会面临较大的偏差风险，在这种情况下采用另一种统计方法更为适宜：即使用MASS包中的rlm函数来进行稳健回归分析
Robust Linear Models（Resilient Linear Models）是一种被广泛认可的线性模型类型，在面对数据集中的异常情况时表现出色。该模型类被设计用于提高在数据不完整或存在噪声的情况下仍能保持良好性能的能力。
mod3<-MASS::rlm(n~wday*term,data=daily)
daily%>%add_residuals(mod3,"resid")%>%ggplot(aes(date,resid))+geom_ref_line(h=0)+geom_line()
六、自然样条法拟合年度平滑曲线模型
> library(splines)
> mod<-MASS::rlm(n~wday*ns(date,5),data=daily)
> daily%>%data_grid(wday,date=seq_range(date,n=13))%>%add_predictions(mod)%>%ggplot(aes(date,pred,color=wday))+geom_line()+geom_point()
通过这样我们也可以看到周六的模式与其他天的模式很不一样。