別說我不教你賣鑽石---基於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倍多，這個後期可能需要優化一下，不過不想那麼多了。所以我覺得我應該拋棄多元線性模型和回歸樹，使用隨機森林模型，所以以後要預測鑽石的價格就使用這個模型；

參考文獻

《數據挖掘與R語言》

《R語言核心技術手冊》

度娘

大家也可以加小編微信：tswenqu，進R語言交流群。