案例 | 基于R语言钻石价格预测

1.1问题描述和目标

因为钻石的价格定价取决于重量，颜色，刀工等影响，价格该如何制定合理，为公司抢占市场制定价格提供依据。

1.2数据说明

这里我使用的是R语言里面数据集diamonds,如果看这本《ggplot2:数据分析与图形艺术》应该对这个数据都不会太陌生。该数据集收集了约54000颗钻石的价格和质量的信息。每条记录由十个变量构成，其中有三个是名义变量，分别描述钻石的切工，颜色和净度；

carat:克拉重量

cut:切工

color:颜色

clarity:净度

depth:深度

table:钻石宽度

以及X，Y,Z

1.3数据加载到R中

由于数据集是R语言自带的，所以我们只要输入下面的命令行查看数据前六行。

head(diamond)

1.3.1构建数据

sample_date <- sample(2,nrow(diamonds),replace=TRUE,prob=c(0.7,0.3))
trai_date <- diamonds[sample_date==1,]#训练数据
test_date <- diamonds[sample_date==2,]#测试数据

1.4数据可视化和概括

因为我们目前没有该问题领域足够多的信息，所以我们先要了解一下这些数据的统计特性是一种比较好的方式，它方便我们更好理解问题。

首先获取概括性统计描述

summary(trai_date)

对于名义变量它给出了每个可能取值的频数，例如,在刀工上ideal等级比其他等级刀工的钻石更多，其他

以下代码回执出关于钻石深度的一个分部

library(car)
par(mfrow=c(1,2))
hist(trai_date$depth,prob=T,xlab=" ")
lines(density(trai_date$depth))
qq.plot(trai_date$depth)

可以看得出来钻石深度在62左右是在最多，分部服从一个类正态分布。Q-Q图上看，数据完全是不服从正态分布的，因为它呈现的是一个曲线，不是一个直线。

这时候我们在使用lattice包绘制箱线图，看看在颜色上的分布是否有区别

library(lattice)
bwplot(depth~color,data=trai_date)

在各个颜色上钻石的深度差异其实不大。

1.5缺失值的处理

数据中出现缺失值的情况其实非常普遍，这会导致在一些不能处理缺失值的分析方法无法应用，一般我们遇到缺失值的数据时候，可以运用几种常见的策略。

~将含有缺失值的案例剔除

~根据变量之间的相关关系填补缺失值

~根据案例之间的相似性填补缺失值

~使用能够处理缺失值数据的工具

这里由于数据集中不存在缺失值，所以以上方法不讲了，请原谅我的坑爹。等下次遇到再讲该如何处理。

1.6各个属性相关性探索

因为我们想要知道各个变量之间的相关性到低是怎么样的，这样子对我们的建模的时候考虑选择模型的时候或者选择变量都有一个很好的参考，这时候我们有个叫cor函数可以计算各个变量之间的相关系数，并产生一个相关系数矩阵

cor(trai_date[,-c(2,3,4)])#这里剔除了三个名义变量

这里为了使其输出结果更加的友好，我们使用symnum函数改善输出结果，这里我们可以看得出钻石的深度貌似和其他变量相关性都不是很强，而钻石重量克拉数却和价格，X,Y,Z相关度特别高，这时候我们对模型的各个变量都有个大致的了解了；是时候展现真正的技术了，那就改考虑模型的选择了

symnum(cor(trai_date[,-c(2,3,4)]))

1.7获取预测模型

因为我们主要是的研究目的是预测，预测测试数的钻石价格；不过从数据结构和数据分布上来看，我们可以使用回归模型和随机森林两类预测模型模型；在回归类的模型中我们可以考虑使用多元线性回归和回归决策树两种模型，到时候我们在建立一个评估模型的函数看哪个模型的预测误差小

1.7.1多元线性回归

这里我们使用Lm函数对数据进行拟合，预测变量是价格，因此我们先初步对多元线性回归模型的一个探索先

lm_model <- lm(price~.,data=trai_date)

summary(lm_model)

得到的模型结果很NICE，调整后的可决系数高达0.9194，这说明了模型可以解释百分九十一的方差，各个变量的显著性也很高，大部分都是三颗星以上的高度显著；先说一下R是怎么处理那三个名义变量的，当像上面一样进行建模的时候，R会生成一组辅助变量，对每一个有K个水平的因子变量，R会生成K-1个辅助变量，辅助比那辆的值为0或者1，当辅助变量的值为1，表示该因子出现，同时表明其他所有辅助变量值为0，以上结果汇总；所以从上图结果我们可以看得出来cut变量生成了cut.l,cut.q,cut.c,cut.4:这4个辅助变量。

我们要提出那些不相关的变量，一个个剔除确实是有些麻烦，这个时候我们选择通过对初始模型用向后剔除法得到一个新的模型

step_lm_model <- step(lm_model)

上图结果展示，模型的AIC值是下降了，也剔除了变量Y,模型是稍微精简了点，这时候我们在看一下模型的整体结果如何；

summary(step_lm_model)

我看到这个结果的时候其实我是拒绝的，因为。。。。；因为模型的结果没变化，尔康我很心塞，苍天不公，好吧，不过因为可决系数很高，所以我们还不能放弃，这时候通过模型诊断看看

先用PLOT函数画出诊断图。

par(mfrow=c(2,2))
plot(step_lm_model)

弱弱解释着这四幅图，解释不好还请各位拍砖头，

左上：代表的残差值和拟合值的拟合图，如果模型的因变量和自变量是线性相关的话，残差值和拟合值是没有任何关系的，他们的分布应该是也是在0左右随机分布，从结果上看，还算是一条直线分布，在后面结尾的时候有几个离群点；

   右上:代表正态QQ图，说白了就是标准化后的残差分布图，如果满足正态假定，那么点应该都在45度的直线上，若不是就违反了正态性假定，开始和结尾是的角度数我不敢恭维，不过我们考虑加个非线性项进去；

   左下：位置尺度图，主要是检验是否同方差的假设，如果是同方差，周围的点应该随机分布，从结果上看，一条直线毕竟无法容纳太多点

   右下：主要是影响点的分析，叫残差与杠杆图，鉴别离群值和高杠杆值和强影响点，说白了就是对模型影响大的点。

加了一个二次项后发现模型结果纹丝不动，意思也就是说么怎么改变；到这里我也该放弃了，毕竟强扭的瓜不甜，强追的女孩受伤重；

1.7.2回归树 

这时候我们需要加载包rpart,然后通过rpart函数构建回归树

library(rpart)
 tree_model <-rpart(price~.,data=trai_date)

在稍微解释在一下这个结果吧，其实已经有写博客介绍过这个结果了，第二行包括了一些信息，包括了节点编号，描述，观察值的数目，偏差和预测值；

对模型进行可视化，这里就不需要我博客课上写过的用maptree包里面的draw.tree函数去画图,这里就还是用一下

library(maptree)
draw.tree(tree_model)

这时候为了防止过度拟合，我们需要对模型进行剪枝，就是偏差的减少小于某一个给定的限定值时候

这里的因为选择不详细介绍，因为篇幅有限，老衲也想早点写完；这时候我们需要确定和计算每个节点的参数值cp,这个参数CP值就是决定函数rpart在构建树的时候如何选择，因此在这里我们生成各个树节点的情况,使用rsq.rpart打印结果

rsq.rpart(tree_model)

这时候我们要选择出合适的CP值来构建复杂性和效果性达到一个平衡点，选择plotcp函数帮助我们选择cp值，根绝下图，选择cp在0.1的时候最佳。

plotcp(tree_model)

这时候我们通过使用prune函数对初始模型进行剪枝，然后得到结果

tree_model2 <- prune(tree_model,cp=0.1)

1.7.3随机森林

写了那么多，还是咬咬牙写一下这个模型，随机森林模型工作原理就是将各个学习模型的结果组合在一起，如果是分类属性，根据投票选出最优结果，如果是连续型，像价格这样的就话就求平均；

选取合适的变量数

因为合适的变量是能够有效的降低误差，因此我们需要写个程序遍历一下，数据集总共有十个变量

m =ncol(trai_date)
for (i in 1:(m-1)){
  test_model <- randomForest(price~.,data=trai_date,mtry=i,importance=TRUE,
                             proximity=TRUE)
  mse <- mean(test_model$mse)
  print(mse)
}

这里模型告诉我们要6个，注意是6个！！因为这时候价格是连续型变量，所以只能要均方残差，如果是字符型变量也就是名义型变量的话就要使用err

选择合适的NTREE值

ntree就是随机森林的的决策树数量，设置过低话预测误差过高，而NTREE过高的话又会提升模型的复杂度，降低效率；初始设置为500，查看结果

forest_model <- randomForest(price~.,data=trai_date,mtry=6,ntree=500)
plot(forest_model)

当子树在100左右的时候误差基本没什么变化，因此我们选择Ntree=100

forest_model_final <- randomForest(price~.,data=trai_date,mtry=6,ntree=100)

这个模型能解读数据的方差占97.81%，看来模型很成功。NICE

1.8模型的评价和选择

在这里我使用一个简单的方法来对模型拟合的判断，我写一个函数，估计每个模型的训练数据和测试数据的均方根误差，当然你也可以自己写的方法，这里不限制

代码如下

 cal_rms_error <- function(model,train,test,yval){
   train.yhat <- predict(object=model,newdata=train)
   test.yhat <- predict(object=model,newdata=test)
   train.y <- with(train,get(yval))
   test.y <- with(test,get(yval))
   train.err <- sqrt(mean((train.yhat-train.y)^2))
   test.err <- sqrt(mean((test.yhat-test.y)^2))
   c(train.err=train.err,test.err)
 }

说明一下这个函数，第一个是模型对象，第二是训练数据，第三是测试数据，第四预测变量，它会输出训练数据和测试数据的均方根误差

首先判断一下多元线性回归模型，看看是不是我们的真爱

cal_rms_error(step_lm_model,trai_date,test_date,'price')

在判断判断一下回归树，突然发现误差挺大的，算了，不考虑它

 cal_rms_error(tree_model2,trai_date,test_date,'price')

最后我们判断I一下随机深林的模型

 cal_rms_error(forest_model_final,trai_date,test_date,'price')

这个结果告诉我，什么才叫做真爱，不过测试数据是训练数据的2倍多，这个后期可能需要优化一下，不过不想那么多了。所以我觉得我应该抛弃多元线性模型和回归树，使用随机森林模型，所以以后要预测钻石的价格就使用这个模型；

【直播课程推荐】

点击“阅读原文”，即可进入学习界面哦~~~