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() #獲取粗房屋名 不過，這就遇到了一個問題，在提取『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 test
D = 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()函式則會返回這些資料集中的一個（預設）或多個。
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)) > 2
GVIF 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，說明不存在多重共線情況。