深圳二手房回归分析

#setwd("/Users/alexandrashaw/Documents/商务数据分析/homework")                    #工作路径

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.4.4
library(gridExtra)


par(family='STKaiti')
dat0=read.csv("mydatashenzhen2.csv",header=T, fileEncoding = 'GBK') #读入清洗过后的数据
dat0=dat0[,-1]                  #去掉第一列序号
n=dim(dat0)[1]                  #样本量
summary(dat0)                   #查看数据基本描述


dat0$price=dat0$danjia/10000            #价格单位转换成万元


#调整城区和楼层的因子水平顺序，以便作图输出美观
dat0$district=factor(dat0$district, levels=c("南山","福田","宝安","龙华","罗湖","龙岗","坪山"))
dat0$floor=factor(dat0$floor,levels=c("高层","中层","低层"))

######对dat0做描述分析

#因变量直方图
par(mfrow=c(1,1))   #画1*2的图
#hist(dat0$danjia, xlab="单位面积房价(万元/平方米)",ylab="频数",main="",col="lightblue")

ggplot(dat0, aes(danjia)) + geom_histogram(fill ="lightblue" ) + xlab("单位面积房价(万元/平方米)") + ylab("频数") + ggtitle("")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
 
#查看房价最高和最低的两条观测
dat0[dat0[,1]==min(dat0[,1]),]

#城区对房价的分组箱线图
#boxplot(danjia~district,data=dat0,col="lightblue",ylab="单位面积房价(万元/平方米)")   

ggplot(dat0, aes(district, danjia)) + geom_boxplot(fill ="lightblue" ) + ylab("单位面积房价(万元/平方米)") 
 
#par(mfrow=c(1,2))  #画1*2的图
#地铁、学区的分组箱线图  #没有地铁附近的房子
#boxplot(danjia~subway,data=dat0,col=c("lightblue","orange"), names=c("非地铁","地铁"),ylab="单位面积房价(万元/平方米)")    
#boxplot(danjia~school,data=dat0,col=c("lightblue","orange"), names=c("非学区","学区"),ylab="单位面积房价(万元/平方米)")    


dat0$schoolname <- ifelse(dat0$school==1, "学区","非学区")
ggplot(dat0, aes(schoolname, danjia, fill = schoolname)) + geom_boxplot() + ylab("单位面积房价(万元/平方米)")   + xlab("") + theme(legend.position =  "none")
 
#par(mfrow=c(2,2))  #画2*2的图
#卧室数、厅数、楼层、朝向的分组箱线图
# boxplot(danjia~room,data=dat0,col="lightblue",ylab="单位面积房价(万元/平方米)",xlab="卧室数")   
# boxplot(danjia~hall,data=dat0,col="lightblue",ylab="单位面积房价(万元/平方米)",xlab="厅数")        
# boxplot(danjia~floor,data=dat0,col="lightblue",names=c("高","中","低"),xlab="楼层",ylab="单位面积房价(万元/平方米)")  
# boxplot(danjia~chaoxiang,data=dat0,col="lightblue",ylab="单位面积房价(万元/平方米)",xlab="朝向")   

p1 <- ggplot(dat0, aes(factor(room), danjia)) + geom_boxplot(fill ="lightblue") + ylab("单位面积房价(万元/平方米)") + xlab("卧室数") 
p2 <- ggplot(dat0, aes(factor(hall), danjia)) + geom_boxplot(fill ="lightblue") + ylab("单位面积房价(万元/平方米)") + xlab("厅数") 
p3 <- ggplot(dat0, aes(factor(floor), danjia)) + geom_boxplot(fill ="lightblue") + ylab("单位面积房价(万元/平方米)") + xlab("楼层") 
p4 <- ggplot(dat0, aes(chaoxiang, danjia)) + geom_boxplot(fill ="lightblue") + ylab("单位面积房价(万元/平方米)") + xlab("朝向") 

grid.arrange(p1,p2,p3,p4, nrow = 2)
 
#par(mfrow=c(1,1))  #画1*1的图
#房屋面积和单位面积房价的散点图
#plot(dat0$area,dat0$danjia,xlab="面积(平方米)",ylab="单位面积房价(万元/平方米)")
ggplot(dat0, aes(area, danjia)) + geom_point() + ylab("单位面积房价(万元/平方米)") + xlab("面积(平方米)") 
 
#par(mfrow=c(2,4))          #画2*4的图
# ###分城区的学区房分组箱线图
# boxplot(danjia~school,data=dat0[dat0$district=="南山",],main="南山")
# boxplot(danjia~school,data=dat0[dat0$district=="福田",],main="福田")
# boxplot(danjia~school,data=dat0[dat0$district=="宝安",],main="宝安")
# boxplot(danjia~school,data=dat0[dat0$district=="龙华",],main="龙华")
# boxplot(danjia~school,data=dat0[dat0$district=="罗湖",],main="罗湖")
# boxplot(danjia~school,data=dat0[dat0$district=="龙岗",],main="龙岗")
# boxplot(danjia~school,data=dat0[dat0$district=="坪山",],main="坪山")


dat <- subset(dat0, district %in% c("南山", "福田","宝安", "龙华","罗湖","龙岗","坪山"))
ggplot(dat, aes(schoolname, danjia, fill = schoolname)) + geom_boxplot() + facet_wrap(~district) +  ylab("单位面积房价(万元/平方米)") +   theme(legend.position =  "none")
 

#导入数据
dat0 <- read.csv("mydatashenzhen2.csv")
dat0=dat0[,-1]					
n=dim(dat0)[1]					
summary(dat0)					

#建立模型
dat0 = dat0[,-ncol(dat0)]

full = lm(danjia ~ ., data = dat0)
summary(full)

#可以看到模型拟合有NA的变量，需要进行变量筛选
#筛选结果
back.aic  = step(full, trace = FALSE)
summary(back.aic )


#模型诊断:

#查看共线性

library(car)
vif(back.aic)

#所有VIFS 值小于5， 没有共线性


#再看回归残差诊断图

par(mfrow = c(2,2))

plot(back.aic, which = c(1, 2, 3, 5))


par(mfrow = c(1,1))

plot(back.aic, which = 4)


#可以从左上的残差图看到，残差呈现方差变大的情况，这意味着残差的方差不等，模型假设不成立，同时COOK距离有明显的强影响点


#下面考虑变换变量

full.log = lm(log(danjia) ~ . - subway, data = dat0)
summary(full.log )

#再查看共线性和残差图

vif(full.log)

#所有VIFS 值小于5， 没有共线性


#再看回归残差诊断图

par(mfrow = c(2,2))

plot(full.log, which = c(1, 2, 3, 5))

par(mfrow = c(1,1))

plot(full.log, which = 4)

#可以看到模型残差方差没有显著变大了，但仍然受到异常点，高杠杆点和影响点 因此去掉它们

drop1 = which(rstandard(full.log) > 2.5)
drop2 = which(hatvalues(full.log) > 2 * mean(hatvalues(full.log)))
drop3 = which(cooks.distance(full.log) > 4 / length(cooks.distance(full.log)))
drop = unique(c(drop1,drop2, drop3))

full.log.refit = lm(log(danjia) ~ . - subway, data = dat0[-drop,])
summary(full.log.refit)
par(mfrow = c(2,2))
plot(full.log.refit, which = c(1, 2, 3, 5))

par(mfrow = c(1,1))

plot(full.log.refit, which = 4)


#此时我们看到残差图左上图（或者左下图） 方差没有增大趋势，方差一致，且图上没有曲线特征，线性性也成立，
#右上图的正态QQ图说明模型残差已经近似正态分布了，点近似组成一条直线了。
#同时我们看到右下角杠杆图没有点显著偏离，且COOK距离图也没有点有明显大的
#COOK距离，因此此时模型的假设都得到满足了 我们可以用此模型进行预测

newdata = data.frame(area = 95, floor = "中层", hall = 3, room = 3, school = 1,
                     chaoxiang = "南向", year = 2018, district = "福田", subway = 0)

pred <- predict(full.log.refit, newdata = newdata, interval = "confidence")
pred <- exp(pred)
pred
#可以看到 一个2018年的95平三房中层位于福田南向的学区房预测为73862.58单价，
#95% 置信区间为[70521.03 77362.47]