乾貨--線性回歸模型與CART樹的比較

作者：劉順祥 公眾號：每天進步一點點2015 (微信ID：lsxxx2011)

配套教程：手把手教你做文本挖掘 https://edu.hellobi.com/course/181

有關CART演算法的理論這裡不再贅述，可參考《淺談C5.0與CART演算法的比較--理論理解》，線性回歸的理論部分也不過多講解，可以參考我之前寫的文章《R語言下的線性回歸模型》和《基於R語言的線性回歸模型診斷》。但我需要強調的是，線性回歸模型的夠著需要滿足幾個前提假設，即：

1）殘差項的期望為0，標準差為sigma的正態分布；

2）殘差項之間是獨立的；

3）殘差項的方差應該滿足齊性；

4）Y變數與X變數之間存在線性關係；

5）X變數之間不存在多重共線性；

之所以需要這些假設，是因為線性回歸的參數確定是通過最小二乘法得到的。如果線性回歸的這些前提假設不能得到保證，那出來的結果也是一種偽回歸，故需要針對這些假設前提給予相應的解決辦法。

線性回歸之5個前提假設的解決辦法分別是：

1）對Y變數作變換，如倒數化、對數化、開根號等，一般可通過BOX-COX變換直接處理；

2）對於非時間序列，一般認為樣本間是獨立的，而針對時間序列的樣本，需要進行杜賓瓦森檢驗，如果數據不滿足獨立性的話，可通過Y變數進行變換，類似於1）中的方法；

3）方差齊性的檢驗可以使用巴特萊特檢驗法，如果不滿足方差齊性，仍然需要對Y變數進行變換，類似於1）中的方法；

4）對於線性關係的檢驗，可以通過皮爾遜檢驗，發現Y變數與各個X變數之間的相關係數，或者通過成分殘差圖來判斷；

5）方差膨脹因子可以鑒別X變數之間是否存在多重共線性，對於高度多重共線性的變數要進行適當的刪除；

實戰

接下來我們就以實戰的方式來比較這兩種用於解決回歸問題的演算法，數據來源於Kaggle網站(https://www.kaggle.com/harlfoxem/housesalesprediction)，該數據集的目的是用來預測房價。首先，我們來了解一下該數據集：

# 載入所需的第三方包

library(GGally)

library(car)

library(rpart)

# 讀取數據，並進行探索性數據分析

house <- read.csv(file.choose())

# 查看數據屬性

dim(house)

names(house)

該數據集一共包含了21613個觀測和21個變數，這些變數包含了房子的價格、卧室數量、洗手間數量、房子的面積、樓層、建築時間、是否重修等。這麼多的變數，到底哪些變數跟房價具有更高的相關性呢？

這裡先刪除一些對Y變數沒有實質意義的自變數，如id、date、zipcode、lat、long。

# 剔除不需要的變數

house.new <- subset(house, select = -c(id,date,zipcode,lat,long))

計算房價與其餘變數的相關係數，由於R自帶的cor函數返回的是所有變數之間的兩兩相關係數，我們這裡不妨自定義一個函數，僅探索Y變數與其餘自變數之間的相關係數：

# 構造自定義函數

mycorr <- function(df, y.index = NULL, y.names = NULL){

if(is.null(y.index) == TRUE)

index = which(names(df) == y.names)

else

index = y.index

# 僅提取出Y變數與自變數的相關係數

test = cor(df)[-index, index]

# 將test結果轉化為矩陣

res = matrix(test, dimnames = list(names(df)[-index],

names(df)[index]))

# 將矩陣轉化為數據框

corr = data.frame(res)

# 給數據框添加變數名，存放各個自變數

corr$variables = row.names(corr)

# 刪除數據框的行名稱

row.names(corr) = NULL

# 將相關係數降序排序

corr = corr[order(corr[,1],decreasing = TRUE),]

# 返回結果

return(corr)

}

mycorr(df = house.new,y.names = price)

發現相關性較高的是sqft_living、grade、sqft_above、sqft_living15和bathrooms五個變數，接下來我們就利用這5個變數來做多元線性回歸分析。

# 篩選相關的變數

house.model <- subset(house.new,

select = c(price,sqft_living,

grade,sqft_above,

sqft_living15,bathrooms))

接下來我們先探索房價的分布情況，看看其是否服從正態分布（一個前提假設，等價於殘差項服從正態分布）。

# 繪製房價的直方圖

hp <- hist(house.new$price,

breaks = (max(house.new$price) - min(house.new$price))/100000,

freq = FALSE,

col = steelblue, xlab = 房價, ylab = 頻率,

main = 房價的直方圖)

lines(density(house.new$price),

col = red, lty = 2, lwd = 2)

x <- seq(min(house.new$price),

max(house.new$price), length = 1000)

y <- dnorm(x, mean = mean(house.new$price),

sd = sd(house.new$price))

lines(x = x, y = y, col = green, lty = 2, lwd = 2)

# 添加圖例

legend(locator(), legend = c(核密度曲線,正態分布曲線),

lty = c(2,2), lwd = c(2,2),

col = c(red,green), bty = n)

很顯然房價並不服從正態分布，後面需要對數據進行BOX-COX變換，以滿足線性回歸的前提假設。

接下來再來看看各個變數之間的詳細情況：

ggpairs(house.model)

這個矩陣圖的對角線反應了某個變數的核密度圖，即實際的分布情況；下三角反應的是某兩個變數間的散點圖；上三角反應的是某兩個變數間的相關係數。

在探索完Y變數和自變數後，我們進一步開始線性模型的搭建，具體可見下方：

# 構建多元線性回歸模型，並對模型的前提假設做檢驗

# 將數據集拆分為訓練集和測試集

set.seed(1234)

index <- sample(1:nrow(house.model),

size = 0.75*nrow(house.model))

train <- house.model[index,]

test <- house.model[-index,]

# 建模

fit1 <- lm(price ~ ., data = train)

summary(fit1)

針對上方的模型反饋結果，我們大致解釋一下：模型中的各個自變數的參數都是通過顯著性檢驗的，同時模型的顯著性檢驗也滿足。但並不代表模型是OK的，因為可能存在偽回歸，這就需要我們對模型的假設前提作一一的檢驗。

# 模型變數間的多重共線性檢驗

vif(fit1)

這個結構能說明什麼呢？根據先驗知識，當0<VIF<10，不存在多重共線性；當10≤VIF<100，存在較強的多重共線性；當VIF≥100，存在嚴重多重共線性。所有，這5個自變數間並不存在較強的多重共線性。

# 殘差的正態性檢驗

# 由於樣本量超過5000，故我們使用K-S檢驗

ks.test(fit1$residuals, pnorm,

mean(fit1$residuals),

sd(fit1$residuals))

上圖也同樣說明了模型的殘差不服從正態性檢驗，那麼該如何處理呢？接下來我們採用BOX-COX檢驗對Y變數進行變換。

powerTransform(fit1)

也許這裡的結果你並不明白，我們先把核對錶擺出來，你大致就可以搞清楚了：

根據變換核對錶，我們發現-0.0043更接近於0，故只需對Y變數作對數變換即可。

# 變換並重新建模

fit2 <- lm(log(price) ~ ., data = train)

summary(fit2)

從上面的結果可知，模型的參數和整個模型也都通過了顯著性檢驗，接下來我們還需要進一步做其他方便的假設檢驗。

# 模型的線性關係檢驗，繪製成分殘差圖

crPlots(fit2)

從圖中可知，變數bathrooms與price之間並不存在明顯的線性關係，而幾乎是一條水平直線，故考慮刪除該變數。

fit3 <- update(fit2, . ~ . - bathrooms)

summary(fit3)

模型結果如下：

# 殘差同方差性檢驗

ncvTest(fit3)

從結果總可知，雖然p值小於0.05，但也比較解決0.05了，說明可以放寬並接受方差齊性的條件。

最後，我們通過一張診斷圖來了解一下模型情況

opar <- par(no.readonly = TRUE)

par(mfrow = c(2,2))

plot(fit3)

par(opar)

左上圖反應了模型的殘差在0附近隨機波動，且不存在某種趨勢，說明殘差的期望為0；

右上圖反應的是模型殘差的QQ圖，即正態性檢驗圖，幾乎所有的點都在虛線上，說明殘差服從正態分布；

左下圖的紅線也沒有明顯的趨勢，說明殘差的方差滿足齊性的假設；

右下圖則反應的數據的異常信息，有個別點屬於異常值，如14號、27號等樣本；

# 預測

pred <- predict(object = fit3, newdata = test[,c(-1,-6)])

使用多元線性回歸模型對房價做了預測之後，我們再來使用CART演算法進行房價的預測：由於CART演算法對數據沒有過多的要求，故我們使用所有的變數對Y變數進行預測。

# 將原始數據進行拆分，構造訓練集和測試集

set.seed(1234)

index <- sample(1:nrow(house.new),

size = 0.75*nrow(house.new))

train2 <- house.new[index,]

test2 <- house.new[-index,]

# 構建CART演算法模型

fit4 <- rpart(price ~ .,data = train2)

# 預測

pred2 <- predict(fit4,newdata = test2[,-1])

模型建完了，預測值也做出來了，如何比較CART演算法和多元線性回歸模型之間的優劣呢？接下來我們使用RMSE(均方根誤差)來衡量。

多元線性回歸模型

# 訓練集上的RMSE

rmse.lm.train <- sqrt(mean(fit3$residuals^2))

# 測試集上的RMSE

rmse.lm.test <- sqrt(mean((log(test[,1])-pred)^2))

rmse.lm.train

rmse.lm.test

通過訓練集和測試集的RMSE比較，並沒有發現多元線性回歸模型模型存在過擬合狀態。

CART演算法

# 訓練集上的RMSE

pred3 <- predict(fit4,newdata = train2[,-1])

rmse.cart.train <- sqrt(mean((log(train2[,1])-log(pred3))^2))

# 測試集上的RMSE

rmse.cart.test <- sqrt(mean((log(test2[,1])-log(pred2))^2))

rmse.cart.train

rmse.cart.test

通過訓練集和測試集的RMSE比較，也沒有發現CART模型存在過擬合狀態。

但通過兩個模型的訓練集RMSE的比較，我們發現多元線性回歸模型比CART演算法要更優秀一些。