數據分析項目實戰：大型商場銷售預測（帶你擠進比賽前100名）

背景

此項目為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 RegularFrequency 522 178 8485 195 4824Proportion 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.003592093highest: 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.9900highest: 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 OUT045Frequency 925 1553 1543 1546 880 1559 1550 1548Proportion 0.065 0.109 0.109 0.109 0.062 0.110 0.109 0.109 Value OUT046 OUT049Frequency 1550 1550Proportion 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 2009Frequency 2439 1553 1550 925 1550 1548 1550 1543 1546Proportion 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 SmallFrequency 1553 4655 3980Proportion 0.152 0.457 0.391------------------------------------------------------------------------Outlet_Location_Type n missing distinct 14204 0 3 Value Tier 1 Tier 2 Tier 3Frequency 3980 4641 5583Proportion 0.280 0.327 0.393------------------------------------------------------------------------Outlet_Type n missing distinct 14204 0 4 Value Grocery Store Supermarket Type1 Supermarket Type2Frequency 1805 9294 1546Proportion 0.127 0.654 0.109 Value Supermarket Type3Frequency 1559Proportion 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.6190highest: 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.csvpred.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.csvpred.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的建模速度和效果值得一試，不用將大把的時間花在模型的建模上，而專註於優化和多次嘗試。

以上就是整個分析和建模的過程，如果你有更好的發現，歡迎跟我交流。

感謝您的閱讀。

推薦閱讀：