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

从训练集中去掉异常点

# 4倍以上的为异常点
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 语言中对于特征变量的处理,缺失值的补充等。随后对比了几种特征选择的方法。 从最后的结果来看,可以看到通过领域知识人工选择的变量已经比较逼近算法选择的最后结果。而后续的几种算法在最后的结果上也没有太大的差别。

参考文章

原创声明,本文系作者授权云+社区-专栏发表,未经许可,不得转载。

如有侵权,请联系 yunjia_community@tencent.com 删除。

编辑于

我来说两句

4 条评论
登录 后参与评论

相关文章

来自专栏技术总结

iOS进阶之CAEmitterLayer

2398
来自专栏QQ会员技术团队的专栏

人人都可以做深度学习应用:入门篇(下)

如果这一轮AI浪潮真的会带来新的一轮科技革命,那么我们相信,它也会遵循类似的发展轨迹,逐步发展和走向普及。如果基于这个理解,或许,我们可以通过积极学习,争取成为...

5.3K2
来自专栏逍遥剑客的游戏开发

Direct3D学习(七):DirectX下天空盒子的实现

1164
来自专栏瓜大三哥

视频压缩编码技术(H.264) 之结构

视频的一场或一帧可用来产生一个编码图像。通常,视频帧可分成两种类型:连续或隔行视频帧。在电视中,为减少大面积闪烁现象,把一帧分成两个隔行的场。显然,这时场内邻行...

601
来自专栏老秦求学

从Iris数据集开始---机器学习入门

代码多来自《Introduction to Machine Learning with Python》. 该文集主要是自己的一个阅读笔记以及一些小思考,小总结...

42510
来自专栏大数据文摘

你真的懂TensorFlow吗?Tensor是神马?为什么还会Flow?

2696
来自专栏深度学习自然语言处理

【论文笔记】中文词向量论文综述(二)

一、Improve Chinese Word Embeddings by Exploiting Internal Structure

583
来自专栏数据科学与人工智能

【算法】TextRank算法为文本生成关键字和摘要

TextRank算法基于PageRank,用于为文本生成关键字和摘要。其论文是:

692
来自专栏AI研习社

《哈利·波特》出版二十周年,教大家用神经网络写咒语!

AI 研习社按:不知道你小时候是否梦想过这样的场景,在对角巷的奥利凡德家挑一把魔杖,十一英寸,冬青木、凤凰羽毛,带着它和一只雪白的猫头鹰奔入国王十字车站,从 9...

3498
来自专栏媒矿工厂

HDR关键技术:HEVC/H.265编码方案

前文我们对HEVC的HDR编码优化技术做了介绍,侧重编码性能的提升。本章主要阐述HEVC中HDR/WCG相关的整体编码方案,包括不同应用场景下的HEVC扩展编码...

990

扫码关注云+社区