数据分析项目实战：大型商场销售预测（带你挤进比赛前100名）
数据分析 | 咨询顾问 | 写作爱好者 | Change Maker
背景
此项目为Analiytics Vidhya的竞赛题目，大型商场销售额预测，该项目提供了从不同城市的10家商店中收集的1559种产品的2013年销售数据。目的是通过2013年的销售数据，建立一个销售额预测模型，预测每个产品在特定商店的销售情况。
写此篇文章的目的是记录我分析和思考的整个过程，并与各位分析师朋友一起交流，帮助自己未来持续改善和提升模型效果。
截止到写此文章为止，我的模型预测效果得分目前是前100名。
竞赛截止到2019年1月1号，如果您在零售行业工作，且在寻找一个优秀的模型或者对此类问题感兴趣，这个项目值得大家参与和探索，相信会有收获。
分析工具为R语言，并且最后通过xgboost提升整个模型效果。
你可以通过这个链接访问到竞赛主页 Big Mart Sales III （报名比赛可见数据集）
提出假设：
目标是预测商店的每个产品的销量，对销售额可能的影响因素如下：
与商店有关的假设:
与商品有关的假设:
数据集
变量描述:
探索数据
#加载包
library(ggplot2)
library(grid)
library(gridExtra)
library(caret)
library(randomForest)
library(xgboost)
library(Hmisc)
library(dplyr)
library(rpart)
library(rpart.plot)
#加载数据
train.data <- read.csv("train.csv",stringsAsFactors = F)
test.data <- read.csv("test.csv",stringsAsFactors = F)
test.data$Item_Outlet_Sales <- NA
#合并数据集，特征工程需要，避免重复执行相同代码
all_bms <- rbind(train.data,test.data)
#描述统计量
describe(all_bms)数据描述信息
all_bms
12 Variables 14204 Observations
------------------------------------------------------------------------
Item_Identifier
n missing distinct
14204 0 1559
lowest : DRA12 DRA24 DRA59 DRB01 DRB13, highest: NCZ30 NCZ41 NCZ42 NCZ53 NCZ54
------------------------------------------------------------------------
Item_Weight
n missing distinct Info Mean Gmd .05 .10
11765 2439 415 1 12.79 5.367 5.940 6.655
.25 .50 .75 .90 .95
8.710 12.600 16.750 19.350 20.250
lowest : 4.555 4.590 4.610 4.615 4.635, highest: 21.000 21.100 21.200 21.250 21.350
------------------------------------------------------------------------
Item_Fat_Content
n missing distinct
14204 0 5
Value LF low fat Low Fat reg Regular
Frequency 522 178 8485 195 4824
Proportion 0.037 0.013 0.597 0.014 0.340
------------------------------------------------------------------------
Item_Visibility
n missing distinct Info Mean Gmd .05 .10
14204 0 13006 1 0.06595 0.05554 0.00000 0.01180
.25 .50 .75 .90 .95
0.02704 0.05402 0.09404 0.13765 0.16330
lowest : 0.000000000 0.003574698 0.003589104 0.003591414 0.003592093
highest: 0.313934693 0.321115010 0.323637245 0.325780807 0.328390948
------------------------------------------------------------------------
Item_Type
n missing distinct
14204 0 16
Baking Goods (1086, 0.076), Breads (416, 0.029), Breakfast (186,
0.013), Canned (1084, 0.076), Dairy (1136, 0.080), Frozen Foods (1426,
0.100), Fruits and Vegetables (2013, 0.142), Hard Drinks (362, 0.025),
Health and Hygiene (858, 0.060), Household (1548, 0.109), Meat (736,
0.052), Others (280, 0.020), Seafood (89, 0.006), Snack Foods (1989,
0.140), Soft Drinks (726, 0.051), Starchy Foods (269, 0.019)
------------------------------------------------------------------------
Item_MRP
n missing distinct Info Mean Gmd .05 .10
14204 0 8052 1 141 71.33 42.88 53.23
.25 .50 .75 .90 .95
94.01 142.25 185.86 231.17 249.89
lowest : 31.2900 31.4900 31.8900 31.9558 31.9900
highest: 266.3226 266.4884 266.5884 266.6884 266.8884
------------------------------------------------------------------------
Outlet_Identifier
n missing distinct
14204 0 10
Value OUT010 OUT013 OUT017 OUT018 OUT019 OUT027 OUT035 OUT045
Frequency 925 1553 1543 1546 880 1559 1550 1548
Proportion 0.065 0.109 0.109 0.109 0.062 0.110 0.109 0.109
Value OUT046 OUT049
Frequency 1550 1550
Proportion 0.109 0.109
------------------------------------------------------------------------
Outlet_Establishment_Year
n missing distinct Info Mean Gmd
14204 0 9 0.986 1998 9.347
Value 1985 1987 1997 1998 1999 2002 2004 2007 2009
Frequency 2439 1553 1550 925 1550 1548 1550 1543 1546
Proportion 0.172 0.109 0.109 0.065 0.109 0.109 0.109 0.109 0.109
------------------------------------------------------------------------
Outlet_Size
n missing distinct
10188 4016 3
Value High Medium Small
Frequency 1553 4655 3980
Proportion 0.152 0.457 0.391
------------------------------------------------------------------------
Outlet_Location_Type
n missing distinct
14204 0 3
Value Tier 1 Tier 2 Tier 3
Frequency 3980 4641 5583
Proportion 0.280 0.327 0.393
------------------------------------------------------------------------
Outlet_Type
n missing distinct
14204 0 4
Value Grocery Store Supermarket Type1 Supermarket Type2
Frequency 1805 9294 1546
Proportion 0.127 0.654 0.109
Value Supermarket Type3
Frequency 1559
Proportion 0.110
------------------------------------------------------------------------
Item_Outlet_Sales
n missing distinct Info Mean Gmd .05 .10
8523 5681 3493 1 2181 1844 188.4 343.6
.25 .50 .75 .90 .95
834.2 1794.3 3101.3 4570.1 5522.8
lowest : 33.2900 33.9558 34.6216 35.2874 36.6190
highest: 10306.5840 10993.6896 11445.1020 12117.5600 13086.9648
------------------------------------------------------------------------- Item_Weight 产品重量存在缺失值，需要处理
- Outlet_Size 商店大小存在缺失值，需要处理
- Item_Outlet_Sales 商品销售额的缺失值符合测试集的数目，不用处理
- Item_Fat_Content 商品低脂变量，有LF，low fat，Low Fat，reg，Regular五个类别，从类别描述看可以分为Low Fat和Regular两个变量，表示低脂和常规
- Item_Visibility 商品的展示百分比数据中存在最低值为0，由于训练集的每个商品数据都有对应的销售额，因此这是异常的值，需要处理
- Item_Type 产品类别归类较多，需要进一步看看能够生成新的归类变量
- Item_MRP 给出了商品的标价，能够计算出对应商品的销量（相对于销售额，销量更易于理解和探索）
- Outlet_Establishment_Year 商店的开店年份19xx年，这样的数据不直观，可以调整为年数
清洗数据
1、处理缺失值：商品重量
由于每个商品的重量应该是相同的，所以用商品的重量均值补插缺失的值
#计算商品重量均值，默认移除缺失值
tmp <- aggregate(Item_Weight~Item_Identifier,data = all_bms,FUN = mean)
#将均值赋值给对应缺失的重量值
for (i in which(is.na(all_bms$Item_Weight))) {
all_bms$Item_Weight[i] <- tmp$Item_Weight[tmp$Item_Identifier == all_bms$Item_Identifier[i]]
}
#检查缺失值，处理后没有缺失值
sum(is.na(all_bms$Item_Weight))2、处理缺失值：商店大小
由于商店大小的值跟我们的前面的假设比较相关，商店的面积可能还与商店的类型等有关联，
#检查商店大小与商店类型的关系
prop.table(table(all_bms$Outlet_Size))
tmp2 <- aggregate(Item_Outlet_Sales~Outlet_Identifier+Outlet_Type+Outlet_Size,data = all_bms,FUN = mean)
tmp2
- 结果看出每个商店的类型，大小与平均销售额的关系，似乎商店的大小确实和设想的一样，与商店的类型有比较强的对应关系
- 可以构建一个决策树来预测和补插商店大小的缺失值
- 个人理解商店大小是比较关键的变量，但因为数据集中的商店数量较少，因此在商店层面的区分效果没有商店类型变量那么好
#构建一个决策树来补插缺失值
#非缺失部分作为训练集，缺失部分作为测试集
fit <- rpart(factor(Outlet_Size)~Outlet_Type,data = all_bms[all_bms$Outlet_Size!="",],method = "class")
pred <- predict(fit,all_bms[all_bms$Outlet_Size=="",],type = "class")
all_bms$Outlet_Size[all_bms$Outlet_Size==""] <- as.vector(pred)
#可以看到每个商店对应的大小都被补充完整了
table(all_bms$Outlet_Size,all_bms$Outlet_Identifier)特征工程与EDA
1、创建一个表示销售量的变量
#计算每个商品的销售量，并将所有小数点大于0的数都向上取整
all_bms$Item_Sales_Vol <- round(all_bms$Item_Outlet_Sales/all_bms$Item_MRP+0.5,0)- 为了贴近实际销量，因此将小数点后大于0的数都向上取整
- 销售量将作为模型的预测变量
2、调整Item_Fat_Content错误的因子水平
#低脂商品因为标注差异导致了5个水平的引子
summary(all_bms$Item_Fat_Content)
#将商品正确归类
all_bms$Item_Fat_Content <- as.character(all_bms$Item_Fat_Content)
all_bms$Item_Fat_Content[all_bms$Item_Fat_Content %in%
c("LF","low fat")] <- "Low Fat"
all_bms$Item_Fat_Content[all_bms$Item_Fat_Content %in%
c("reg")] <- "Regular"
table(all_bms$Item_Fat_Content)- 将商品的低脂调整为：Low Fat，Regular两个水平
3、对Item_Type进一步归类，创建新的商品分类变量
#商品的类别可以进一步归类
summary(all_bms$Item_Type)
#通过查看商品ID，可以看到ID开头似乎是商品的大类，DR饮品,DF食物,NC费消耗品，DR和DF为日常消耗品
#ID后的A-Z是商品编码，似乎看不到品牌信息
all_bms$Item_Attribute <- factor(substr(all_bms$Item_Identifier,1,2))
table(all_bms$Item_Attribute)
#由于食物才有低脂的情况，因此对于非消耗品则需要不同区分
all_bms$Item_Fat_Content[all_bms$Item_Attribute == "NC"] <- "Non-Food"
table(all_bms$Item_Fat_Content)- 新建变量Item_Attribute：DR饮品,DF食物,NC费消耗品三个类别
- 由于存在非消耗品，对Item_Fat_Content新增一个因子，Non-Food 非食品
4、将Item_Visibility中为0的商品调整为每个商店中的平均值
#采用每店的平均商品展示百分比来补插为0的异常值
tmp3 <- aggregate(Item_Visibility~Outlet_Identifier,data = all_bms,FUN = mean)
#将商品展示百分比均值数据赋值给为0的值
for (i in which(all_bms$Item_Visibility==0)) {
all_bms$Item_Visibility[i] <- tmp3$Item_Visibility[tmp3$Outlet_Identifier == all_bms$Outlet_Identifier[i]]
}
sum(all_bms$Item_Visibility==0)5、将Outlet_Establishment_Year转化为年数，创建新变量
#2013年收集的数据
all_bms$Outlet_Years <- 2013-all_bms$Outlet_Establishment_Year- 新变量Outlet_Years保存每个商店的成立年数，为离散型变量，打算采用因子类型
6、将所有分类变量转化为因子，消除不存在的因子
cols <- c("Item_Fat_Content","Item_Type","Outlet_Location_Type","Outlet_Type","Outlet_Years","Item_Attribute","Outlet_Identifier")
for (i in cols) {
all_bms[,i] <- factor(all_bms[,i])
}可视化探索各个变量
经过特征工程，初步完成了对我所关心变量的调整，在这一步探索各个变量他们对销售量的表现
1、商店层面的4个因素：Outlet_Type，Outlet_Location_Type，Outlet_Years，Outlet_Size
可视化的结果表明：
- 商店类型因素能够较明显区分每个商店的销量水平
- 商店的地理位置似乎很难看出区别，但其主要原因在于商店数量较少，难以在这个维度区分开来，不能否定该因素的价值
- 商店的年份看，在同个类型的商店，似乎越新的商店有着更好的表现，似乎消费者会更喜欢新开张的商店，而不是因为年代久知名度高的原因，与我们前面的假设有点出入。
- 从开店的趋势看，似乎管理公司更倾向于商品销量表现平稳的类型1商店
- 从商店的大小看似乎也很难看出商品销量的区别，其原因也在于商店数量较少，该变量还是一个比较关键的，因为商店大小跟开店成本相关，在其他预测问题上表现应该很不错，但在此处，可能表现会一般
2、商品层面的因素：Item_Visibility，Item_Weight，Item_MRP，Item_Attribute，Item_Fat_Content，Item_Type（控制商店类型）
可视化结果表明：
- 商品展示百分比，商品价格，商品品类对销量有一定的影响
- 但在商品层面因素，销量的表现差异不如商店层面因素那么明显
- 其他的因素似乎看不到明显的关系
因此，数据建模中，打算采用如下7个变量作为预测变量，进行销量的预测：Outlet_Type，Outlet_Location_Type，Outlet_Years，Outlet_Size，Item_Visibility，Item_MRP，Item_Type
分割数据集
train <- all_bms[!is.na(all_bms$Item_Outlet_Sales),]
test <- all_bms[is.na(all_bms$Item_Outlet_Sales),]
#设置种子，创建训练数据的训练子集和测试子集
set.seed(1234)
ind <- createDataPartition(train$Item_Sales_Vol,p=.7,list = FALSE)
train_val <- train[ind,]
test_val <- train[-ind,]初步建模过程
1、建模准备
#需要建模的变量
myformula <- Item_Sales_Vol ~ Outlet_Type + Item_Visibility + Outlet_Location_Type + Item_MRP + Item_Type + Outlet_Years + Outlet_Size
#创建模型评估RMSE函数
model.rmse <- function(pred,act){
sqrt(sum((act - pred)^2)/
length(act))
}2、决策树（回归）
下面省略了利用增加其他变量进行预测的过程，
#决策树建模
fit.tr <- rpart(myformula,
data = train_val,
method = "anova")
summary(fit.tr)
rpart.plot(fit.tr)
pred <- predict(fit.tr,test_val)
model.rmse(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales)
#创建用于上传评分的测试结果 dtree.csv
pred.test <- predict(fit.tr,test)
submit <- data.frame(Item_Identifier = test.data$Item_Identifier, Outlet_Identifier = test.data$Outlet_Identifier,Item_Outlet_Sales=pred.test*test$Item_MRP)
write.csv(submit, file = "dtree.csv", row.names = FALSE)模型得分：1158.65757134
3、随机森林
set.seed(2345)
fit.rf <- randomForest(myformula,
data = train_val,
ntree=500)
summary(fit.rf)
pred <- predict(fit.rf,test_val)
model.rmse(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales)
#创建用于上传评分的测试结果 rf.csv
pred.test <- predict(fit.rf,test)
submit <- data.frame(Item_Identifier = test.data$Item_Identifier, Outlet_Identifier = test.data$Outlet_Identifier,Item_Outlet_Sales=pred.test*test$Item_MRP)
write.csv(submit, file = "rf.csv", row.names = FALSE)模型得分：1153.04001353
随机森林建模效果有一个明显的提升
4、GBM建模
#控制器,5折交叉验证
Ctrl <- trainControl(method="repeatedcv",number=5,
repeats=5)
set.seed(3456)
fit.gbm <- train(myformula,
data = train_val,
trControl=Ctrl,
method="gbm",
verbose=FALSE
)
summary(fit.gbm)
pred <- predict(fit.gbm,test_val)
model.rmse(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales)
pred.test <- predict(fit.gbm,test)
submit <- data.frame(Item_Identifier = test$Item_Identifier, Outlet_Identifier = test$Outlet_Identifier,Item_Outlet_Sales=pred.test*test$Item_MRP)
write.csv(submit, file = "gbm5cv.csv", row.names = FALSE)模型得评分1152.66471101
gbm效果比随机森林好，当然可以尝试其他参数调整模型，下面的表格可以看到模型中各个变量的重要程度
当然，这样的分数无法让我们进入前100名，而且随机森林和gbm的建模时间有点久，不利用调优。
xgboost建模并进行调优
1、构建稀疏矩阵
由于xgboost模型要求所有的变量都为数值型，因此存在分类变量则需要将分类变量转化为0，1格式的稀疏矩阵
#构建稀疏矩阵
mymatrix <- function(train){
matrix_num <- train[,c("Item_Visibility","Item_MRP")]
matrix_num <- cbind(matrix_num,
model.matrix(~Outlet_Type-1,train),
model.matrix(~Outlet_Location_Type-1,train),
model.matrix(~Outlet_Size-1,train),
model.matrix(~Item_Type-1,train),
model.matrix(~Outlet_Years-1,train)
)
return(data.matrix(matrix_num))
}
#获取每个数据集的稀疏矩阵
xgb.train_val <- mymatrix(train_val)
xgb.test_val <- mymatrix(test_val)
xgb.test <- mymatrix(test)
#生成用于xgboots模型的DMatrix
#注意，预测变量集和响应变量是要分开的
dtrain_val <- xgb.DMatrix(data = xgb.train_val,label=train_val$Item_Sales_Vol)
dtest_val <- xgb.DMatrix(data = xgb.test_val,label=test_val$Item_Sales_Vol)
dtest_sub <- xgb.DMatrix(data = xgb.test)2、初步xgboost建模
#初步xgboost建模
model <- xgboost(data = dtrain_val,nround = 5)
summary(model)
pred <- predict(model,dtest_val)
model.rmse(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales)
xgb.importance(colnames(xgb.train_val),model)
pred.test <- predict(model,dtest_sub)
submit <- data.frame(Item_Identifier = test$Item_Identifier, Outlet_Identifier = test$Outlet_Identifier,Item_Outlet_Sales=pred.test*test$Item_MRP)
write.csv(submit, file = "xgb.csv", row.names = FALSE)评分只有1206.37915379，这比GBM效果要差，因为我们初始设置了循环5次，xgb通过不断的循环来提升模型效果。
3、调优模型
上述结果明显是模型拟合得不够好，因此增加循环次数为10，但为了防止过度拟合，将默认的树模型最大深度6调整为5
model_tuned <- xgboost(data = dtrain_val,
nround = 10,
max.depth = 5
)
summary(model_tuned)
pred <- predict(model_tuned,dtest_val)
model.rmse(pred*test_val$Item_MRP,test_val$Item_Outlet_Sales)
xgb.importance(colnames(xgb.train_val),model_tuned)
pred.test <- predict(model_tuned,dtest_sub)
submit <- data.frame(Item_Identifier = test$Item_Identifier, Outlet_Identifier = test$Outlet_Identifier,Item_Outlet_Sales=pred.test*test$Item_MRP)
write.csv(submit, file = "xgbn10d5.csv", row.names = FALSE)评分：1145.50824952669 模型效果大幅度提升，挤进前100名没问题了
4、微调模型
我尝试到最好，基于目前特征工程后创建的模型，在模型调整循环11次，最大深度为4，可以达到1142.66240303726，处于目前TOP50的状态。
5、查看模型中各个变量的权重
xgb.importance(colnames(xgb.train_val),model_tuned)
结论
- 整个模型也比较容易理解，商店的类型占据绝大部分影响因素，其次就是商店内的商品标价和展示，商品的品类也有一定的影响
- 与建模开始前预测的差不多，商店的地域和商店的大小在此模型中没有突出的贡献，但不能否定这两个变量在同类问题中的重要性，只是因为其对应的商店数据太少
- 我也尝试过直接用销售额进行预测，但效果没有销量的好，这也说明了，为模型提供一个更为简单且易于理解的响应变量对模型的预测效果有一个明显的提升
就建模工具而已，xgboost的建模速度和效果值得一试，不用将大把的时间花在模型的建模上，而专注于优化和多次尝试。
以上就是整个分析和建模的过程，如果你有更好的发现，欢迎跟我交流。
感谢您的阅读。
编辑于 2018-03-27 22:06
