Kaggle实战：House Prices: Advanced Regression Techniques（下篇）
原创Kaggle实战：House Prices: Advanced Regression Techniques（下篇）
原创
作者：郭小发
接上一篇《Kaggle 实战-House Prices: Advanced Regression Techniques（上篇）》
初步模型
这次题目给的自变量有很多，我们需要从中挑选对房价影响最大的变量。 我们的思路是先人工挑选一些对房价影响比较重要的因素，然后再慢慢的添加新的变量来看是否会改变模型的精度。
以国内的房价为例，影响房价的因素主要是房子面积、房子所在的区域、小区等，房龄、房型（小高层、多层、别墅等）、特殊场景（地铁房、学区房等）、装修等也会影响价格。 这个数据是美国的房屋信息，不过基本的影响因素应该差不多。
我们先来简单的撸出一个模型来，选择如下变量：
- LotArea 房子的面积
- Neighborhood 城市街区 用来初步代替 区域、小区
- Condition1 Condition2 附近的交通情况
- BldgType 房屋类型 独栋别墅、联排别墅
- HouseStyle 房子的层数
- YearBuilt 房子建造的年份
- YearRemodAdd： 房子的改造年份
- OverallQual： 房子整体质量，考量材料和完成度
- OverallCond：房子整体条件
变量观察
我们先看一下这些变量和房价之间的关系。
为了方便画图，我们写了两个函数用来画图
# 加载库
library(ggplot2)
# 将对于因子变量画图
plot2_factor <- function(var_name){
source('D:/RData/comm/multiplot.r')
plots <- list()
plots[[1]] <- ggplot(train, aes_string(x = var_name, fill = var_name) ) +
geom_bar() +
guides(fill = FALSE) +
ggtitle(paste("count of ", var_name)) +
theme(axis.text.x = element_text(angle = 90, hjust =1))
plots[[2]] <- ggplot(train, aes_string(x = var_name, y = "SalePrice", fill = var_name) ) +
geom_boxplot() +
guides(fill = FALSE) +
ggtitle(paste( var_name, " vs SalePrice")) +
theme(axis.text.x = element_text(angle = 90, hjust =1))
multiplot(plotlist = plots, cols = 2)
}
# 对于连续数字变量画图
plot2_number <- function(var_name){
source('D:/RData/comm/multiplot.r')
plots <- list()
plots[[1]] <- ggplot(train, aes_string(x = var_name) ) +
geom_histogram() +
ggtitle(paste("count of ", var_name))
plots[[2]] <- ggplot(train, aes_string(x = var_name, y = "SalePrice") ) +
geom_point() +
ggtitle(paste( var_name, " vs SalePrice"))
multiplot(plotlist = plots, cols = 2)
}
# 街区和房价的关系
plot2_factor("Neighborhood")
plot2_number("YearBuilt")
plot2_number("OverallQual")
上图是 Neighborhood 和 SalePrice 之间的对比图。
通过图可以看到不同街区的房价分布还是有很大不同的，这个变量应该很有潜力。
上图是 YearBuilt 和 SalePrice 之间的对比图。
通过图可以看到建造时间越近的房子价格越高。
上图是 OverallQual 和 SalePrice 之间的对比图。
通过图可以看到装修越好的房子价格越高。
查看各变量之间的相关性
# 相关系数画图
library(corrgram)
sel <- c("LotArea","Neighborhood","BldgType","HouseStyle","YearBuilt","YearRemodAdd","OverallQual","OverallCond","MSZoning")
corrgram(train[,sel], order=TRUE, lower.panel=panel.shade, upper.panel=panel.pie, text.panel=panel.txt)
模型训练
先用我们挑选的变量来建立一个 lm 模型，作为我们的 base 模型
# 通过人工选择的变量来构造一个公式
fm.base <- SalePrice ~ LotArea + Neighborhood + BldgType + HouseStyle + YearBuilt + YearRemodAdd + OverallQual + OverallCond
# 训练模型
lm.base <- lm(fm.base, train)
# 查看模型概要
summary(lm.base)
结果如下：
Call:
lm(formula = fm.base, data = train)
Residuals:
Min 1Q Median 3Q Max
-208970 -20882 -2917 15544 351199
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.455e+06 1.850e+05 -7.862 7.42e-15 ***
LotArea 1.084e+00 1.156e-01 9.375 < 2e-16 ***
NeighborhoodBlueste -1.068e+03 2.953e+04 -0.036 0.971141
NeighborhoodBrDale -1.440e+04 1.518e+04 -0.949 0.342806
NeighborhoodBrkSide -1.876e+04 1.278e+04 -1.468 0.142460
NeighborhoodClearCr -2.352e+03 1.332e+04 -0.177 0.859842
NeighborhoodCollgCr -2.917e+04 1.086e+04 -2.685 0.007335 **
NeighborhoodCrawfor 1.747e+04 1.246e+04 1.402 0.161225
NeighborhoodEdwards -2.813e+04 1.165e+04 -2.414 0.015924 *
NeighborhoodGilbert -4.030e+04 1.157e+04 -3.484 0.000508 ***
NeighborhoodIDOTRR -3.357e+04 1.343e+04 -2.499 0.012570 *
NeighborhoodMeadowV 1.338e+04 1.446e+04 0.925 0.354867
NeighborhoodMitchel -2.819e+04 1.196e+04 -2.356 0.018617 *
NeighborhoodNAmes -2.202e+04 1.130e+04 -1.950 0.051426 .
NeighborhoodNoRidge 6.105e+04 1.226e+04 4.980 7.13e-07 ***
NeighborhoodNPkVill 6.340e+03 1.650e+04 0.384 0.700928
NeighborhoodNridgHt 4.876e+04 1.104e+04 4.417 1.08e-05 ***
NeighborhoodNWAmes -2.126e+04 1.166e+04 -1.823 0.068457 .
NeighborhoodOldTown -2.915e+04 1.243e+04 -2.344 0.019194 *
NeighborhoodSawyer -2.575e+04 1.188e+04 -2.168 0.030350 *
NeighborhoodSawyerW -2.224e+04 1.154e+04 -1.927 0.054234 .
NeighborhoodSomerst -1.228e+04 1.093e+04 -1.123 0.261764
NeighborhoodStoneBr 5.984e+04 1.249e+04 4.790 1.84e-06 ***
NeighborhoodSWISU -2.365e+04 1.433e+04 -1.651 0.099024 .
NeighborhoodTimber -1.326e+04 1.236e+04 -1.073 0.283489
NeighborhoodVeenker 2.303e+04 1.555e+04 1.481 0.138905
BldgType2fmCon 1.230e+03 7.413e+03 0.166 0.868218
BldgTypeDuplex -7.231e+02 5.831e+03 -0.124 0.901330
BldgTypeTwnhs -6.675e+04 7.811e+03 -8.546 < 2e-16 ***
BldgTypeTwnhsE -4.916e+04 4.892e+03 -10.049 < 2e-16 ***
HouseStyle1.5Unf -2.835e+04 1.102e+04 -2.573 0.010184 *
HouseStyle1Story -3.981e+03 3.977e+03 -1.001 0.316972
HouseStyle2.5Fin 5.328e+04 1.472e+04 3.619 0.000306 ***
HouseStyle2.5Unf -4.606e+03 1.250e+04 -0.368 0.712613
HouseStyle2Story 4.069e+03 4.205e+03 0.968 0.333393
HouseStyleSFoyer -1.173e+04 7.791e+03 -1.505 0.132424
HouseStyleSLvl -6.438e+03 6.197e+03 -1.039 0.299077
YearBuilt 4.285e+02 8.428e+01 5.084 4.20e-07 ***
YearRemodAdd 3.114e+02 7.505e+01 4.149 3.53e-05 ***
OverallQual 2.849e+04 1.187e+03 24.010 < 2e-16 ***
OverallCond 1.613e+03 1.150e+03 1.402 0.161035
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 38880 on 1419 degrees of freedom
Multiple R-squared: 0.7671, Adjusted R-squared: 0.7605
F-statistic: 116.8 on 40 and 1419 DF, p-value: < 2.2e-16
结果解读
针对模型 summary 之后的结果，我们简单解读一下输出结果含义
残差统计量
Residuals:
Min 1Q Median 3Q Max
-208970 -20882 -2917 15544 351199
线性回归的计算基于一些假设，其中一个假设就是 误差符合相互独立、均值为 0 的正态分布。
从本例可以看出这个残差的中位数为负数，数据整体左偏。其中的 1Q 和 3Q 是第一四分位（first quartile）和第三四分位（third quartile）。残差的最大值和最小值附近对应的记录则可能是异常值。
由于残差代表预测值和真实值之间的差别，也就是说最大值 351199 表示我们预测的最大误差有 35 万美元之多。
仅仅从残差的五数概括上看不出什么关键信息，后续可以通过残差图来检查残差是否符合正态分布的趋势。
回归系数
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.455e+06 1.850e+05 -7.862 7.42e-15 ***
LotArea 1.084e+00 1.156e-01 9.375 < 2e-16 ***
NeighborhoodBlueste -1.068e+03 2.953e+04 -0.036 0.971141
NeighborhoodBrDale -1.440e+04 1.518e+04 -0.949 0.342806
NeighborhoodBrkSide -1.876e+04 1.278e+04 -1.468 0.142460
...
...
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
线性回归拟合完成后得出的回归系数并不是准确的值，而是对于真实回归系数的估计值。
既然是估计值则必然存在误差，上述结果中的
- Estimate 表示回归系数的估计
- Std. Error 表示回归系数的标准误差
- t value 表示假设此回归系数为 0 时的 T 检验值
- Pr(>|t|) 则是上述假设成立的置信度 p-value
P-value 越小则说明假设（假设回归系数为 0）越不容易出现，反过来就是此变量的回归系数不为 0 的几率越大，故此变量在整个回归拟合中作用越显著。一般用置信度 0.05 作为判断依据。
- 最后的三颗星表示此变量显著，星号越多越显著，最多三个。
- 最后一行 Signif. codes 标识着显著标识编码 当 P-value 小于 0.001 时三颗星，小于 0.01 时两颗星，大于 0.05 则认为不太显著。
R2 和调整 R2
Multiple R-squared: 0.7671, Adjusted R-squared: 0.7605
R-squared（判定系数，coefficient of determination） 也称为模型拟合的确定系数，取值 0~1 之间，越接近 1，表明模型的因变量对响应变量 y 的解释能力越强。 Adjusted R-squared 当自变量个数增加时，尽管有的自变量与 y 的线性关系不显著，R square 也会增大。Adjusted R square 增加了对变量增多的惩罚，故我们以 Adjusted R square 为判断模型好坏的基本标准。
本例中 Adjusted R-squared: 0.7605 表示响应变量有 76%的方差被此模型解释了。
模型整体的 F 检验
F-statistic: 116.8 on 40 and 1419 DF, p-value: < 2.2e-16 F 统计量用来检验模型是否显著 假设模型所有的回归系数均为 0，即该模型是不显著的。对此假设做 F 检验，在 p-value 的置信度下拒绝了此假设，则模型为显著的。
在本例中 p-value: < 2.2e-16，远远低于 0.05，所以模型是显著的。
变量选择 - 人工筛选
模型 lm.base 的 Adjusted R-squared:是 0.7605。 从第一个模型的结果看到变量 OverallCond 并不显著，所以我们去掉变量 OverallCond，重新进行拟合。 拟合结果 Adjusted R-squared: 0.7603 和之前相差不大， 并且所有的变量都显著。 故我们将第一个模型定为：
# 初步决定的 lm.base 模型的变量
fm.base <- SalePrice ~ LotArea + Neighborhood + BldgType + HouseStyle + YearBuilt + YearRemodAdd + OverallQual
# 训练模型
lm.base <- lm(fm.base, train)
模型出来了，我们先计算答案提交一次，看看结果如何。
# 用 lm.base 模型预测
lm.pred <- predict(lm.base, test)
# 写出结果文件
res <- data.frame(Id = test$Id, SalePrice = lm.pred)
write.csv(res, file = "D:/RData/House/res_base.csv", row.names = FALSE)
Your submission scored 0.20943
评分标准是测试集上的 root mean squared logarithmic error ，其值越小表示误差越小。 结果惨不忍睹，排名已经 2000 开外了。
初步优化
建立模型之后，我们接下来要诊断一下模型是否符合线性回归的一些假设。 直接通过 plot 给出模型的线性回归诊断图：
# 快速打印残差图、QQ 图等
layout(matrix(1:4,2,2))
plot(lm.base)
- 残差-拟合图（Residuals vs Fitted）
线性回归的计算基于一些假设，其中假设误差是相互独立、均值为 0 正态分布。如果因变量与自变量线性相关的，那么残差的分布应该是正态分布。 通过上图可以看出，残差整体是随机分布在均线 0 值附近的。残差比较大的点很大几率是异常点，需要去除掉。
- 尺度-位置图（Scale-Location Graph）
因变量的方差不随自变量的水平不同而变化，称为同方差性（残差方差不变）。 如果残差满足不变方差假设，那么在尺度-位置图中，水平线周围的点应该随机分布，那么红线应该是比较直的一条线。
- 正态 Q-Q 图（Normal Q-Q）
正态 Q-Q 图是在正态分布对应的值下标准化残差的概率图。如果残差满足相互独立、均值为 0 正态分布的假设，那么图上的点应该落在呈 45 度角的直线上。通过图上可以看到异常值的残差偏离 45 度线比较多。
- 残差与杠杆图（Residuals vs Leverage）
这个图形主要用来鉴别出离群点、高杠杆值点和强影响点。拟合比较好的模型中所有的点都不应该超过 0.5 倍 Cook 距离，也即是不超过图中 0.5 的那根红色点线。
我们通过诊断图看到整体的模型里面有很多的离群点或者异常值，这些异常值会影响模型的整体拟合质量。所以我们下一步则通过 Cook 距离来去除掉所有的异常点。
一般对于 4 倍以上的 Cook 平均距离的点作为异常点的监测标准。
# 通过 cook 距离来查看异常点
cooksd <- cooks.distance(lm.base)
# 画图
plot(cooksd, pch=".", cex=1, main="Influential Obs by Cooks distance") # plot cook's distance
abline(h = 4*mean(cooksd, na.rm=T), col="red") # add cutoff line
text(x=1:length(cooksd)+1, y=cooksd, labels=ifelse(cooksd>20*mean(cooksd, na.rm=T),names(cooksd),""), col="red") # add labels
从训练集中去掉异常点
# ４倍以上的为异常点
influential <- as.numeric(names(cooksd)[(cooksd > 4*mean(cooksd, na.rm=T))])
train <- train[ -influential, ]
我们再查看一下数字变量的分布，其中 SalePrice 做了对数处理后，分布更加趋于正态分布。
# 查看 SalePrice 的分布
layout(matrix(1:2,1,2))
hist(train$SalePrice)
hist(log(train$SalePrice))
因此我们将数字变量都做对数处理后，生成如下的公式
# 新的函数
fm.base <- log(SalePrice) ~ log(LotArea) + Neighborhood + BldgType + HouseStyle + YearBuilt + YearRemodAdd + OverallQual + OverallCond
# 训练模型
lm.base <- lm(fm.base, train)
重新训练后，得到新的的 Adjusted R-squared: 0.8203，相比之前的有所提高。
我们重新提交一下结果，最后的预测结果是 0.17168，误差也相对变小。
至此我们将这个模型作为我们的基础模型。
不同的变量选择方法对比
回归分析中最重要的莫过于变量选择。我们在前面根据我们自己对题目的理解来人工选择了自变量。
下面我们通过一些算法来帮助我们来选择自变量。
这些方法都已经有比较成熟的 R 包来实现。
逐步回归
逐步回归中，模型会一次添加或者删除一个变量，直到达到某个判停准则为止。
向前逐步回归（forward stepwise）每次添加一个预测变量到模型中，直到添加变量不会使模型有所改进为止。
向后逐步回归（backward stepwise）从模型包含所有预测变量开始，一次删除一个变量直到会降低模型质量为止。
我们这里选择向前逐步回归。
#####################################################
# 取空函数和全函数
null=lm(log(SalePrice)~1, data=train)
full=lm(log(SalePrice)~ .-Id , data=train)
# 向前计算
set.seed(999)
lm.for <- step(null, scope=list(lower=null, upper=full), direction="forward")
summary(lm.for)
最后选择的变量如下：
# 最后选择的变量
log(SalePrice) ~ OverallQual + Neighborhood + GrLivArea + BsmtFinSF1 +
OverallCond + YearBuilt + TotalBsmtSF + GarageCars + MSZoning +
BldgType + Functional + LotArea + SaleCondition + CentralAir +
Condition1 + BsmtFullBath + Exterior1st + Fireplaces + YearRemodAdd +
Heating + ScreenPorch + WoodDeckSF + LotFrontage + Foundation +
KitchenQual + BsmtExposure + HeatingQC + SaleType + GarageCond +
KitchenAbvGr + EnclosedPorch + BsmtFinSF2 + X3SsnPorch +
HalfBath + FullBath + LotConfig + Street + GarageArea + OpenPorchSF
LASSO
一般的线性回归使用最小二乘 OLS 进行回归计算很容易造成过拟合，噪声得到了过分的关注，训练数据的微小差异可能带来巨大的模型差异。 而 Lasso 方法使用 L1 正则，解出的参数常常具有稀疏的特征，即很多特征对应的参数会为零，也就淘汰了一些自变量对于因变量的影响。
我们用 R 包 glmnet 来实现 LASSO 算法。
# 安装
install.packages("glmnet")
library(glmnet)
# 准备数据
formula <- as.formula( log(SalePrice)~ .-Id )
# model.matrix 会自动将分类变量变成哑变量
x <- model.matrix(formula, train)
y <- log(train$SalePrice)
#执行 lasso
set.seed(999)
lm.lasso <- cv.glmnet(x, y, alpha=1)
# 画图
plot(lm.lasso)
# 得到各变量的系数
coef(lm.lasso, s = "lambda.min")
#由于 SalePrice 为 NA 无法数组化
test$SalePrice <- 1
test_x <- model.matrix(formula, test)
# 预测、输出结果
lm.pred <- predict(lm.lasso, newx = test_x, s = "lambda.min")
res <- data.frame(Id = test$Id, SalePrice = exp(lm.pred))
write.csv(res, file = "D:/RData/House/res_lasso.csv", row.names = FALSE)
随机森林
随机森林能够处理很高维度的数据，并且不用做特征选择，会自动的选取变量，所以直接将所有变量都放进去。
#加载随机森林包
library(randomForest)
library(caret)
#设定种子
set.seed(223)
# 设定控制参数
# method = "cv" -- k 折交叉验证
# number -- K 折交叉验证中的 K， number=10 则是 10 折交叉验证
# repeats -- 交叉验证的次数
# verboseIter -- 打印训练日志
ctrl <- trainControl(method = "cv", number = 10, repeats = 20, verboseIter = TRUE)
# 训练模型
lm.rf <- train(log(SalePrice)~ .-Id, data = train, method = "rf", trControl = ctrl, tuneLength = 3)
# 输出结果
write_res(lm.rf, test, 'rf')
# 输出结果
lm.pred <- predict(lm.rf, test)
res <- data.frame(Id = test$Id, SalePrice = exp(lm.pred))
write.csv(res, file = "D:/RData/House/res_rf.csv", row.names = FALSE)
梯度下降 GBDT
GBDT 全称为 Gradient Boosting Decision Tree,它是一种基于决策树(decision tree) 实现的分类回归算法。它和随机森林一样都是模型组合的一种，都是将简单的模型组合起来，效果比单个更复杂的模型好。 组合的方式不同导致算法不同，随机森林用了随机化方法，而 GBDT 则使用了 Gradient Boosting 的方法。
我们用 R 包 gbm 来实现 GBDT 算法。
# 安装包
install.packages("gbm")
# 训练模型
lm.gbm <- train(log(SalePrice)~ .-Id, data = train, method = "gbm", trControl = ctrl)
# 输出结果
lm.pred <- predict(lm.gbm, test)
res <- data.frame(Id = test$Id, SalePrice = exp(lm.pred))
write.csv(res, file = "D:/RData/House/res_gbm.csv", row.names = FALSE)
汇总结果
我们上面使用了不同的算法来对特征选择，我们提交答案的最后结果如下：
结论
这篇文章主要根据实例演示了 R 语言中对于特征变量的处理，缺失值的补充等。随后对比了几种特征选择的方法。 从最后的结果来看，可以看到通过领域知识人工选择的变量已经比较逼近算法选择的最后结果。而后续的几种算法在最后的结果上也没有太大的差别。
参考文章
- Regression Analysis by Example . Samprit Chatterjee, Ali S. Hadi
- R in Action, Second Edition . Robert I. Kabacoff
- Machine Learning With R Cookbook . Yu-Wei, Chiu (David Chiu)
- Glmnet Vignette： http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html#qs
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。