線性回歸和時間序列分析北京房價影響因素可視化案例|附代碼數(shù)據(jù)
全文鏈接：http://tecdat.cn/?p=21467
最近我們被客戶要求撰寫關(guān)于北京房價的研究報告，包括一些圖形和統(tǒng)計輸出。
在本文中，房價有關(guān)的數(shù)據(jù)可能反映了中國近年來的變化
目的
人們得到更多的資源（薪水），期望有更好的房子
人口眾多
獨生子女政策：如何影響房子的幾何結(jié)構(gòu)？更多的臥室，更多的空間
我核心的想法是預(yù)測房價。然而，我不打算使用任何arima模型；相反，我將使用數(shù)據(jù)的特性逐年擬合回歸。
結(jié)構(gòu)如下：
數(shù)據(jù)準(zhǔn)備：將數(shù)值特征轉(zhuǎn)換為分類；缺失值
EDA:對于數(shù)值特征和分類特征:平均價格與這些特征的表現(xiàn)
建模：
分割訓(xùn)練/測試給定年份的數(shù)據(jù)：例如，在2000年分割數(shù)據(jù)；根據(jù)這些數(shù)據(jù)訓(xùn)練回歸模型
然后，在2016年之前的所有新年里，預(yù)測每套房子的價值。
用于驗證的度量將是房屋的平均價格（即每年從測試樣本中獲得平均價格和預(yù)測值）
數(shù)據(jù)準(zhǔn)備
我們對特征有了非常完整的描述：
url：獲取數(shù)據(jù)（字符）的urlid:id（字符）Lng：和Lat坐標(biāo)，使用BD09協(xié)議。（數(shù)字）Cid：社區(qū)id（數(shù)字）交易時間：交易時間（字符）DOM：市場活躍日。（數(shù)字）關(guān)注者：交易后的人數(shù)。（數(shù)字）總價：（數(shù)值）價格：按平方計算的平均價格（數(shù)值）面積：房屋的平方（數(shù)字）起居室``數(shù)（字符）客廳``數(shù)（字符）廚房：廚房數(shù)量（數(shù)字）浴室數(shù)量（字符）房子高度建筑類型：包括塔樓（1）、平房（2）、板塔組合（3）、板（4）（數(shù)值）施工時間裝修：包括其他（1）、粗（2）、簡單（3）、精裝（4）（數(shù)值）建筑結(jié)構(gòu)：包括未清（1）、混合（2）、磚和木（3）、磚混凝土（4）、鋼（5）和鋼-混凝土復(fù)合材料（6）（數(shù)值）梯梯比：同層居民數(shù)與電梯數(shù)量的比例。電梯有（1）或沒有電梯（0）（數(shù)值）五年期：業(yè)主擁有不到5年的財產(chǎn)（數(shù)字）
數(shù)據(jù)清理、特征創(chuàng)建
從最初的數(shù)據(jù)看：
從網(wǎng)址上，我發(fā)現(xiàn)它有位置信息，如chengjiao/101084782030。同樣，一個簡單的regexp進(jìn)行省特征提取。
另一個大的數(shù)據(jù)準(zhǔn)備工作是轉(zhuǎn)換一些數(shù)字特征，比如地鐵，地鐵站附近的房子編碼為1，相反的情況編碼為0。
還有很大一部分DOM缺失。我既不能在建模中使用這個特性，也不能刪除NA，但它也會減小數(shù)據(jù)幀的大小。
#從網(wǎng)址中提取省份
?sapply(df$url, function(x) strsplit(x,'/')[[1]][4])
檢查缺失
#缺失數(shù)據(jù)圖
?ggplot(data = .,aes(x = V2, y = V1)) + geom_tile(aes(fill = value )) +
?
如上所述，DOM的很大一部分丟失了。我決定先保留這個特性，然后用中間值來填充缺失的值（分布是非常傾斜的）
否則，buildingType和communityAverage（pop.）中只有幾個缺少的值，我決定簡單地刪除這些值。事實上，它們只占了約30行，而整個數(shù)據(jù)集的數(shù)據(jù)量為300k+，因此損失不會太大。
下面我簡單地刪除了我以后不打算使用的特征。
ifelse(is.na(df$DOM),median(df$DOM,na.rm=T),df$DOM)
點擊標(biāo)題查閱往期內(nèi)容
數(shù)據(jù)分享|R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機(jī)森林算法預(yù)測心臟病
左右滑動查看更多
01
02
03
04
用于將數(shù)字轉(zhuǎn)換為類別的自定義函數(shù)
對于某些特征，需要一個函數(shù)來處理多個標(biāo)簽，對于其他一些特征（客廳、客廳和浴室），轉(zhuǎn)換非常簡單。
df2$livingRoom?<-?as.numeric(df2$livingRoom)
似乎buildingType具有錯誤的編碼數(shù)字值：
buildingTypecount0.04840.12530.25020.33350.37510.42910.500150.66711.000845412.0001373.000597154.000172405NaN2021
由于錯誤的編碼值和NA的數(shù)量很少，因此我將再次丟棄這些行
df2$renovationCondition <- sapply(df2$renovationCondition, ionCondition)
df2$buildingStructure <- sapply(df2$buildingStructure, makeStructure)
df2$elevator <- ifelse(df2$elevator==1,'has_elevator','no_elevator')
缺失值檢察
# 缺失數(shù)據(jù)圖
df2 %>% is.na %>% melt %>%
?ggplot(data = .,aes(x = Var2, y = Var1)) + geom_tile(aes(fill = value)) + ?scale_fill_manual(values = c("grey20","white")) + theme_minimal(14) +
kable(df?%>% group_by(constructionTime) %>% summarise(count=n()) %>% arrange(-count) %>% head(5))constructionTimecount200421145200319409NA19283200518924200614854
?
df3?<-?data.frame(df2?%>%?na.omit())
插補(bǔ)后的最終檢查
any(is.na(df3))## [1] FALSE
探索性分析
由于有數(shù)字和分類特征，我將使用的EDA技術(shù)有：
數(shù)值：相關(guān)矩陣
分類：箱線圖和地圖
我們必須關(guān)注價格（單位價格/單位價格）以及總價格（百萬元）
totalPrice將是回歸模型的目標(biāo)變量。
數(shù)值特征
corrplot(cor(
?df3 ?, ?tl.col='black')
評論
totalPrice與communityAverage有很強(qiáng)的正相關(guān)關(guān)系，即人口密集區(qū)的房價較高
totalPrice與客廳、衛(wèi)浴室數(shù)量有一定的正相關(guān)關(guān)系。
至于面積變量，我們看到它與上述變量也有很強(qiáng)的相關(guān)性：這是有道理的，因為如果房子的面積大，可以建造更多的房間（顯而易見）。
其他一些有趣的相關(guān)性：communityAverage與建筑時間呈負(fù)相關(guān)，這意味著在人口密集區(qū)建房所需的時間更短
分類特征
地圖
中國三級（省）地圖
我看了看城郊，它位于北京附近，所以我過濾了那個特定省份的地圖
ggplot() +
?geom_polygon(data = shapefile_test,aes(x = long, y = lat, group = group),
BeijingLoc <- data.frame('Long'=116.4075,'Lat' = 39.904)
建筑結(jié)構(gòu)
makeEDA('buildingStructure'?)
磚木結(jié)構(gòu)的房屋是最昂貴的，幾乎是其他類型房屋的兩倍
點擊標(biāo)題查閱往期內(nèi)容
R語言用線性回歸模型預(yù)測空氣質(zhì)量臭氧數(shù)據(jù)
左右滑動查看更多
01
02
03
04
建筑類型
makeEDA('buildingType'?)
平房是最昂貴的
裝修條件
電梯
價格對電梯的依賴性非常小
住宅的分布與這一特征是相對相等的。
地鐵
價格對地鐵站附近的依賴性非常小。
住宅的分布與這一特征是相對相等的。
是否滿_五年_
makeFeatureCatEDA('fiveYearsProperty',?length(unique(df3$fiveYearsProperty)))
對于是否擁有不到5年房產(chǎn)來說，價格的依賴性確實很小
就這一特征而言，房子的分布是相對平等的
區(qū)域
回歸模型
策略
從tradeTime中提取年份和月份
按年度和月份分組，得到房屋的數(shù)量和均價拆分?jǐn)?shù)據(jù)集：
對于年[2010-2017]=在這組年上訓(xùn)練并運行回歸模型
對于>2017年：逐月對測試樣本并預(yù)測平均價格
平均價格總覽
首先我們需要看看我們想要預(yù)測什么
df3$year <- year(df3$tradeTimeTs)
df3$month <- month(df3$tradeTimeTs)df3 %>% filter(year>2009) %>% group_by(monthlyTrad) %>%
?summarise(count=n(), mean = mean(price)) %>%
?ggplot(aes(x=monthlyTradeTS, y= mean)) +
平均價格上漲至2017年中期，然后迅速下降
同時，房屋數(shù)量隨著價格的上漲而增加，而且現(xiàn)在房屋交易的數(shù)量也隨著價格的上漲而減少。
準(zhǔn)備訓(xùn)練/測試樣本
我在2017-01-01拆分?jǐn)?shù)據(jù)。對于所有樣本，我需要把分類特征變成偽變量。
df_train <- data.frame(df ?%>% filter(year>2009 & year<2017))
df_test <- data.frame(df %>% filter(year>=2017))
as.data.frame(cbind(
?df_train %>% select_if(is.numeric) %>% select(-Lng, -Lat, -year, -month), ?'bldgType'= dummy.code(df_train$buildingType), ?'bldgStruc'= dummy.code(df_train$buildingStructure), ?'renovation'= dummy.code(df_train$renovationCondition), ?'hasElevator'= dummy.code(df_train$elevator),
在這一步中，我只訓(xùn)練一個線性模型
regressors<-c('lm')
Control <- trainControl(method = "cv",number = 5, repeats=3)for(r in regressors){
? ?cnt<-cnt+1
? ? res[[cnt]]<-train(totalPrice ~., data = train ,method=r,trControl = ?Control)
?r^2在0.88左右，不錯。讓我們看看細(xì)節(jié)。
訓(xùn)練精度
g1<-ggplot(data=PRED,aes(x=Prediction,y=True)) + geom_jitter() + geom_smooth(method='lm',size=.5) +
? ?#計算指標(biāo)
? ?mse <- mean((PRED$True-PRED$Prediction)^2)
? ?rmse<-mse^0.5
? ?SSE = sum((PRED$Pred - PR
## [1]?"MSE: 15952.845934 RMSE : 126.304576 R2 :0.795874"
所以看起來殘差還不錯（分布是正態(tài)的，以0為中心），但對于低價格來說似乎失敗了。
?訓(xùn)練和測試樣本的預(yù)測與時間的關(guān)系
基本上與上述相同，但我將重復(fù)預(yù)測所有月份的訓(xùn)練數(shù)據(jù)
我的目標(biāo)指標(biāo)是平均房價。
訓(xùn)練是在10多年的訓(xùn)練樣本中完成的，因此逐月查看預(yù)測將非常有趣。
# 訓(xùn)練樣本->訓(xùn)練精度
for (i in 1:length(dates_train)){
? ? current_df <- prepareDF(current_df)
? ? current_pred <- mean(predict(res[[1]],current_df))
#運行測試樣本-->測試精度
for (i in 1:length(dates_test)){
? ? current_df <- prepareDF(current_df)
? ?current_pred <- mean(predict(res[[1]],current_df)) RES %>% reshape2::melt(id=c('date','split')) %>%
?ggplot(aes(x=date,y=value)) + geom_line(aes(color=variable, lty=split),size=1) +
預(yù)測對于2012年之后的數(shù)據(jù)確實非常好，這可能與有足夠數(shù)據(jù)的月份相對應(yīng)
改進(jìn)
地理位置作為特征
下面是一個有趣的圖；它顯示了每個位置的總價格。在二維分布的中心，價格更高。
這個想法是計算每個房子到中心的距離，并關(guān)聯(lián)一個等級/分?jǐn)?shù)
BeijingLoc <- data.frame('Long'=116.4075,'Lat' = 39.904)
df3 %>% ggplot(aes(x=Lng,y=Lat)) + geom_point(aes(color=price),size=.1,alpha=.5) ?+
?theme(legend.position = 'bottom') +
本文摘選?《?R語言線性回歸和時間序列分析北京房價影響因素可視化案例?》?，點擊“閱讀原文”獲取全文完整資料。
點擊標(biāo)題查閱往期內(nèi)容
向量自回歸（VAR）模型分析消費者價格指數(shù) (CPI) 和失業(yè)率時間序列
Matlab用BUGS馬爾可夫區(qū)制轉(zhuǎn)換Markov switching隨機(jī)波動率模型、序列蒙特卡羅SMC、M H采樣分析時間序列
Matlab創(chuàng)建向量自回歸（VAR）模型分析消費者價格指數(shù) (CPI) 和失業(yè)率時間序列
Stata廣義矩量法GMM面板向量自回歸 VAR模型選擇、估計、Granger因果檢驗分析投資、收入和消費數(shù)據(jù)R語言時變向量自回歸（TV-VAR）模型分析時間序列和可視化
R語言用向量自回歸（VAR）進(jìn)行經(jīng)濟(jì)數(shù)據(jù)脈沖響應(yīng)研究分析
R語言arima，向量自回歸（VAR），周期自回歸(PAR)模型分析溫度時間序列
R語言VAR模型的不同類型的脈沖響應(yīng)分析
R語言隨機(jī)搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
R語言時變參數(shù)VAR隨機(jī)模型
R語言估計時變VAR模型時間序列的實證研究分析案例
R語言向量自回歸模型(VAR)及其實現(xiàn)
R語言實現(xiàn)向量自回歸VAR模型
R語言估計時變VAR模型時間序列的實證研究分析案例
Python和R用EWMA，ARIMA模型預(yù)測時間序列
R語言用LASSO，adaptive LASSO預(yù)測通貨膨脹時間序列
Python中的ARIMA模型、SARIMA模型和SARIMAX模型對時間序列預(yù)測
R語言arima，向量自回歸（VAR），周期自回歸(PAR)模型分析溫度時間序列
【視頻】Python和R語言使用指數(shù)加權(quán)平均(EWMA)，ARIMA自回歸移動平均模型預(yù)測時間序列
