R的爬虫和回归模型案例-以北京自如房租价格为例

一、背景

爬虫向来不是R的专长，但并不代表R在此方面一事无成。正好在学习R的rvest爬虫包，不如边学边做，自己来做一个案例。
作为一名消费者，自如不错的管理服务和靠谱的房源使得在帝都想省事儿的我们，即使花高于周边其他竞争者的价格，也愿意租住自如。不过，只要是住过自如一两年或看过自如几次房子的小伙伴应该会注意到，自如的房价与我们实际所看到房间的自己预估价格，并不是都契合的，总是会存在被高估和被低估的房子。那么，作为一名消费者，在与商家的博弈中，消费者手上有更为透明的信息无疑会有更强的议价能力。
所以，本案例的目的在于，通过分析现有房租价格和房间信息之间的关系，建立基于这些信息的定价体系回归模型。

二、数据收集

自如作为一个现在最大的房屋租赁中介，已经在其官网上放置了海量的数据，所以，本案例分析的所有数据自然全部来自北京自如官网。
在开始介绍爬虫之前，先就自如官网的数据信息的几个特点说明一下：

a、该网站信息是处于实时更新状态，本次爬虫获得的数据是2017年1月5日14:26获得的；
b、该网站的租房信息也不是那一时刻自如所有的信息，自如有还多待租但还没有上网的房源（是的，我的根据是每次找房时候自如管家总会神秘的对我说，『除了网上的，我们还有几个没上网的房间，要不看一下』）；
c、本次分析的目标信息是自如合租月付这类房子，整租、日付、自如寓等以后有精力再说；
d、本次分析的数据只包括北京城六区（东城、西城、朝阳、海淀、丰台、石景山）外加昌平的自如合租月付数据，一方面是个人找房方位的喜好，一方面是其网页设置的原因，待会儿会提到。

1、首先是找数据源。

2、其次是不断观察总结各个数据节点，找到目标数据所在节点。

3、再次是对照包的参数，形成相应的爬虫代码。

library(rvest)
library(stringr)
library(XML)
library(xml2)WebSpider <- function(m){url <- str_c(cp,"?p=",m)web <- read_html(url,encoding = "UTF-8")#抓取网页信息name_rough <- web %>% html_nodes("h3")  %>%  html_text() #获取粗房屋名area_rough <- web %>% html_nodes("h4")  %>% html_text() #提取区位price_rough <- web %>% html_nodes("p.price")  %>% html_text() #提取价格price <- str_extract(price_rough, "[0-9]+") %>% as.numeric()#提取精确价格detail <- web %>% html_nodes("div.detail")  %>%  html_text() #提取其他信息#合并成数据框data.frame(name_rough,area_rough,forward,mate_num,location,price,detail)
}

不过，这就遇到了一个问题，在提取『detail』一块时，无法再往下有区分的抓取了，因此造成最后抓到的数据是这样的：

只能在Excel里操作整理完成了。
4、接着是观察翻页规律，然后遇到了一个坑。原以为之后的页码不过是http://www.ziroom.com/z/nl/z3.html的基础上加上/1、/2……..，我依照这个思路抓了一番，也获得了将近4000多条数据，原以为这大概就是全部吧。但我仔细看数据才发现，这样下来的基本都是房山、大兴和通州等的数据，基本没有城六区的，城六区的只有选了区域选项后才会出现：

dc <- "http://www.ziroom.com/z/nl/z3-d23008614.html"
xc <- "http://www.ziroom.com/z/nl/z3-d23008626.html"
cy <- "http://www.ziroom.com/z/nl/z3-d23008613.html"
hd <- "http://www.ziroom.com/z/nl/z3-d23008618.html"
ft <- "http://www.ziroom.com/z/nl/z3-d23008617.html"
sjs <- "http://www.ziroom.com/z/nl/z3-d23008623.html"
cp <- "http://www.ziroom.com/z/nl/z3-d23008611.html"

这样一来，只有逐区的来进行翻页爬了。为此，只能选定部分区域来做分析了。

results_cp <- data.frame()
for(m in 1:118){ #118为昌平区信息的总页码 results_cp <- rbind(results_cp,WebSpider(m))#合并单个区每一次循环输出的数据
}
#依次重复获得7个区的数据
results <- rbind(results_cp, results_cy,results_dc,results_ft,results_hd,results_sjs,results_xc) #将所有各区数据的合并

最后，因为#3中遇到的问题，不得不在excel里先进行一轮数据清洗和整理，大致就是筛选、分列和去重的工作，得到这样的数据格式，最终获得12250条数据（去重前14060）并导入R。

三、数据清洗和考察

在在excel里清理过一次后，还是在R里围绕分析目标做一些数据整理。

library(readxl)
library(stringr)
library(dplyr)
library(psych)
library(ggplot2)
library(nortest)
library(gridExtra)
library(mice)
library(VIM)
library(corrplot)
library(DMwR)
library(car)zr_data <- read_excel("~/Desktop/Über R/Excesise/ziroom-excel-V4.xlsx")zr_data$method_of_pay <- str_replace(zr_data$method_of_pay, "[M]", "1")
zr_data$method_of_pay <- str_replace(zr_data$method_of_pay, "[D]", "2")%>% as.numeric()
zr_mp <- subset(zr_data, method_of_pay == "1")#subset the data of pay monthly
zr_dp <- zr_data[zr_data$method_of_pay == "2", ]#subset the data of pay dayly
zr_mp_c <- filter(zr_mp, method_of_rent == "合")#subset the data of cotenant from monthly paying

最终筛选获得符合本次目标的数据集，自如合租月付（而非日付）的9712条数据：zr_mp_c。通过summary大致了解数据后，为数据集里的数据做一个归类。

#group each variable by the data types
norminal_data <- list(zr_mp_c$forward, zr_mp_c$district, zr_mp_c$area, zr_mp_c$layout)
ordinal_data <- list(zr_mp_c$stock)
ratio_data <- list(zr_mp_c$price, zr_mp_c$area_space, zr_mp_c$distance_from_ss, zr_mp_c$stock)

单变量探索

正态性检验

首先是房租价格、房屋面积、地铁距离远近和所在楼层数几个连续性变量的正态性检验。

h1 <- ggplot(zr_mp_c, aes(x = price, y = ..density..)) + geom_histogram(binwidth = 100) + geom_density()+labs(title = "房租价格")+theme(text = element_text(family = "STSong"))
h2 <- ggplot(zr_mp_c, aes(x = area_space, y = ..density..)) + geom_histogram(binwidth = 1) + geom_density()+labs(title = "房屋面积")+theme(text = element_text(family = "STSong"))
h3 <- ggplot(zr_mp_c, aes(x = distance_from_ss, y = ..density..)) + geom_histogram(binwidth = 50) + geom_density()+labs(title = "距离地铁远近")+theme(text = element_text(family = "STSong"))
h4 <- ggplot(zr_mp_c, aes(x = stock, y = ..density..)) + geom_histogram(binwidth = 1) + geom_density()+labs(title = "楼层分布")+theme(text = element_text(family = "STSong"))grid.arrange(h1, h2, h3, h4, nrow = 2)

以上图片大致可见，房租价格和房屋面积虽然不是标准的正态分布，但符合正态性的特征可大致确定，距离地铁远近和楼层分布就不太好确定。
接下来是用函数来做更为确切的检验。我之前一般会用shapiro和ks.test来检测，不过这次做这个案例也学到了一个点：ks.test可用于大样本但要求样本中不能出现相同值；而shapiro检验又不适合大样本的检验（3~5000之间）。为此，我找到了nortest包的lillie.test，其是ks.test的修正。
代码和结果如下：

sapply(ratio_data, lillie.test)> sapply(ratio_data, lillie.test)[price]
statistic 0.05587243
p.value   2.151243e-80
method    "Lilliefors (Kolmogorov-Smirnov) normality test"                                    [area_space]
statistic 0.1050441
p.value   8.539925e-294
method    "Lilliefors (Kolmogorov-Smirnov) normality test"                                      [distance_from_ss]
statistic 0.09777697
p.value   3.632065e-238
method    "Lilliefors (Kolmogorov-Smirnov) normality test"  lillie.test(zr_mp_c$stock)> lillie.test(zr_mp_c$stock)Lilliefors (Kolmogorov-Smirnov) normality testD = 0.19381, p-value < 2.2e-16

结合函数的正态性检验，基本可以确定四个变量都是符合正态分布的。不过结合观察这几个结果可以看到，距离地铁远近和楼层分布并无异常。房租价格和房屋面积在右侧有明显长尾，说明可能有异常数据。
调出数据来看，这几个不管是价格或面积高得异常的数据，回到原网站查看，却是存在，排除整理分析过程错误。至于是不是数据录入错误，根据对于北京区位的了解，也觉得可能性很小。最后作为可信数据留下了。

缺失值考察

接着是缺失值的考察。在开始之前，先需要挑选最终进入模型构建的变量，并对其中的离散型变量进行因子化处理。

 #construct as the final data frame
zr_mp_c_a <- subset(zr_mp_c, select = c("price", "area_space", "distance_from_ss", "stock", "forward", "district", "area", "layout"))
#form the new stock variable
table(stock2 <- cut(zr_mp_c_a$stock,breaks = c(1,6,12,20,33), labels = c("low","low midlle","high middle","high")))low  low midlle high middle        high 4697        2136        1508         630
zr_mp_c_a$stock2 <- stock2#transform the factor variables
zr_mp_c_a$forward <- factor(zr_mp_c_a$forward)
zr_mp_c_a$district <- factor(zr_mp_c_a$district)
zr_mp_c_a$area <- factor(zr_mp_c_a$area)
zr_mp_c_a$layout  <- factor(zr_mp_c_a$layout)

主要变量的缺失值状况：

> sum(is.na(zr_mp_c_a$price))
[1] 0
> sum(is.na(zr_mp_c_a$distance_from_ss))
[1] 598
> sum(is.na(zr_mp_c_a$stock))
[1] 6
> sum(is.na(zr_mp_c_a$stock2))
[1] 6
> sum(is.na(zr_mp_c_a$area_space))
[1] 13
> sum(is.na(zr_mp_c_a$forward))
[1] 5
> sum(is.na(zr_mp_c_a$district))
[1] 0
> sum(is.na(zr_mp_c_a$area))
[1] 1
> sum(is.na(zr_mp_c_a$layout))
[1] 0

可以看到在诸多指标中，距离地铁远近这个指标的缺失值最多，占总值的6%左右，其余的缺失值都很少。接着可以借助VIM包的matrixplot看一下缺失值的分布状况。

matrixplot(zr_mp_c_a)

通过分布可以看到，距离地铁远近这个指标的缺失值分布较为集中，因为数据整体是以区为单位排布的，基本可以确定是某几个区（大致是昌平和海淀吧）在这个指标的信息录入上有问题。这也就意味着有必要处理缺失值，不然会影响到特定指标的特征。
本次原本打算用knnImputation方便快捷地来插补的，不过一直报invalid “time” arguments的错误，遂继续尝试用更为智能的mice。其大致原理是

mice包提供了多种先进的缺失值处理方法。它使用一种不同寻常的方法来进行两步插值：首先利用mice函数建模再用complete函数生成完整数据。mice(df)会返回df的多个完整副本，每个副本都对缺失的数据插补了不同的值。complete()函数则会返回这些数据集中的一个（默认）或多个。

关于mice的缺失值插值处理介绍，可移步在R中填充缺失数据-mice包和R语言-缺失值处理。

micemod <- mice(zr_mp_c_a[1:4]))
zr_mp_c_a[1:4] <- complete(micemod)

在完成插补后，也就获得了一个较为完全的数据集。

变量间关系的探索

关于变量间相关性的探索，可以使用psych包的corr.test一次性的完成相关程度和显著性的检验。

corr.test(zr_mp_c_a[ ,1:4])

可以看到，几个连续性变量之间除了楼层数量未通过显著性检验完外，其他三个都具有显著性。同时，各变量间的相关程度也不算高，没有明显的变量共线情况。

四、模型建构

模型对比

准备工作做的差不多了，就直接开始上模型。

fit1 <-  lm(price ~ ., data = zr_mp_c_a)
#compare the influence of "area" and "districk"
fit2 <-  lm(price ~ area_space + distance_from_ss + forward + stock + area + layout, data = zr_mp_c_a)
fit3 <-  lm(price ~ area_space + distance_from_ss + forward + stock +district + layout, data = zr_mp_c_a)
#drop "stock" and "area"
fit4 <-  lm(price ~ area_space + distance_from_ss + forward + district+ layout, data = zr_mp_c_a)

fit1放入了所有预测变量，虽然获得0.7086的R方值，但「area」变量庞杂的变量类别使得模型散失了可解读性，所以在fit2和fit3中，对比了这两个变量剔除后的变化。此外，在双变量对比时发现「stock」与「price」的相关性没有显著性，在fit4里将其剔除。最终，fit4的模型特征如下：

模型评估

绘制fit4 的诊断图

opar <- par(no.readonly = TRUE)
par(mfrow = c(2,2))
plot(fit4)
par(opar)

对照模型拟合情况诊断图，逐一来看：

a、首先是正态性（右上图）大部分都落到了直线上，可大致认为分布通过检验；
b、独立性：诊断图无法获得，那就借助car包里的durbinWatsonTest函数来判断

   > durbinWatsonTest(fit4)  lag Autocorrelation D-W Statistic p-value    1   0.07321289996    1.85339171       0  Alternative hypothesis: rho!= 0

p值不显著说明无自相关性，独立性可获得认可；

c、线性：（左上图）的基本算是一条直线，可算是直线关系，不必做进一步的处理。
d、位置尺度图：（左下图）存在较多离散的点，可以考虑对响应变量进行对数处理。

fit4 <- lm(log(price) ~ area_space + distance_from_ss + forward + district+layout, data = zr_mp_c_a)

再次对比诊断图：

响应变量的离散情况有所改善。

诊断多重共线情况

> vif(fit4)GVIF Df GVIF^(1/(2*Df))
area_space       1.136749561  1     1.066184581
distance_from_ss 1.126070012  1     1.061164461
forward          1.198440060  9     1.010107448
district         1.259109540  6     1.019385909
layout           1.172449041 17     1.004690222
> sqrt(vif(fit4)) > 2GVIF    Df GVIF^(1/(2*Df))
area_space       FALSE FALSE           FALSE
distance_from_ss FALSE FALSE           FALSE
forward          FALSE  TRUE           FALSE
district         FALSE  TRUE           FALSE
layout           FALSE  TRUE           FALSE

总体来说几个预测变量的方差膨胀因子都 <2，说明不存在多重共线情况。