干货--线性回归模型与CART树的比较

有关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算法要更优秀一些。

今天我们就先分享到这里，如有任何问题都可以给我留言，谢谢。下一期我们将道道分类算法中的knn来龙去脉。

每天进步一点点2015

学习与分享，取长补短，关注小号！