R语言波士顿房价分析

最近帮GF用R语言做了一个波士顿房价的多元线性回归的作业，因为也是初次学习，可能内容会有一些问题，希望大家一起学习交流。
R用的是3.5.3版本，下载都是来自China.Lanzhou.
##########程序如下所示###########
###安装更新

install.packages(“installr”)
library(installr)
updateR()
#####安装包###
install.packages(“readr”)
install.packages(“funr”)
install.packages(“openxlsx”)
install.packages(“dplyr”)
install.packages(“caret”)
install.packages(“readxl”)
install.packages(“psych”)
install.packages(“plyr”)
install.packages(“ggplot2”)
install.packages(“rio”)
install.packages(“abind”)

##########加载包#####

library(readr)
library(funr)
library(openxlsx)
library(dplyr)
library(caret)
library(readxl)
library(psych)
library(plyr)
library(ggplot2)
library(zoo)
library(car)
############ 导入数据 ##########
path <- get_script_path()# 获取脚本路径
path = “C:/Users/shaon/Desktop/LWY/大作业”
setwd(path)
model_data <- read_csv(“Boston-house-prices.csv”)
str(model_data)#可以看到导入的数据
summary(model_data$MEDV)

hist(model_data$MEDV)#可以看到目标数据的总体分布

############ 数据分析#######
###形成相关系数###

cor(model_data[c(“CRIM”,“ZN”,“INDUS”,“CHAS”,“NOX”,“RM”,“AGE”,“DIS”,“RAD”,“TAX”,“PTRATIO”,“BRATIO”,“LSTAT”,“MEDV”)])
###散点图矩阵###
pairs(model_data[c(“CRIM”,“ZN”,“INDUS”,“CHAS”,“NOX”,“RM”,“AGE”,“DIS”,“RAD”,“TAX”,“PTRATIO”,“BRATIO”,“LSTAT”,“MEDV”)])

###改进后的散点图矩阵###
pairs.panels(model_data[c(“CRIM”,“ZN”,“INDUS”,“CHAS”,“NOX”,“RM”,“AGE”,“DIS”,“RAD”,“TAX”,“PTRATIO”,“BRATIO”,“LSTAT”,“MEDV”)])

#####基于数据训练模型,得到目标函数的系数#####
lmrelationship<-lm(MEDV~CRIM+ZN+INDUS+CHAS+NOX+RM+AGE+DIS+RAD+TAX+PTRATIO+BRATIO+LSTAT,data=model_data)
summary(lmrelationship)

可以看到R^2只有0.7406和0.7338，比较小，同时变量中有很多与目标变量相关度较低的，这里我们删去一些相关度不高的参数，只留下NOX,RM,DIS,RAD,PTRATIO,LSTAT，然后重复上述工作。

lmrelationship2<-lm(MEDV~NOX+RM+DIS+RAD+PTRATIO+LSTAT,data=model_data)
summary(lmrelationship2)

发现R^2的变化不大，说明删除的变量还算合理，这里如果希望删除更多变量或者提高R可以从不相关度从大到小逐个删除来比较。

#############异方差的检验#############
###G-Q检验###

res.test<-residuals(lmrelationship2)
library(lmtest)
gqtest(lmrelationship2)###p值很小时拒绝H0，认为上诉公式有相关性，存在异方差

###BPtest检验###

bptest(lmrelationship2)####p值很小时认为存在异方差

################修正异方差##################
异方差的修正关键在于取什么样的权重，我看网上其他资料采用的权重有1/e,1/e^2,1/x,1/sqrt(x),x ^2等等，这里以权重为1/e为例,

lmrelationship3<-resid(lmrelationship2)^2###将精简后模型的残差赋给lmrelationship3
lmrelationship4<-lm(MEDV~NOX+RM+DIS+RAD+PTRATIO+LSTAT,weights=1/abs(lmrelationship3),data=model_data)
summary(lmrelationship4)
deviance(lmrelationship4)###用来求残差平方和RSS

我们可以看到，R^2为0.9994有了很大的改善，RSS=501.2323这个有点大，但还不知道如何解决。
整理公式：MEDV=40.514-19.857NOX+4.0434RM-1.1900DIS+0.0515RAD-1.1300PTRATIO-0.5934LSTAT

###多重共线性检验

library(graphics)
library(car)
vif(lmrelationship4)#VIF小于10则共线性较小。

我们看到各个变量VIF都比较小，说明多重共线性不大。如果需要消除多重共线性可以采用逐步回归，或者手动删除变量。

lmrelationship2<-step(lmrelationship)