数学建模社区-数学中国

标题: 【高级数理统计R语言学习】2 多元线性回归 [打印本页]

作者: 1047521767 时间: 2021-10-29 11:44
标题: 【高级数理统计R语言学习】2 多元线性回归
【高级数理统计R语言学习】2 多元线性回归

一、背景: t$ V) ~# W% M0 z" `
数据集展示了X市外来人口的相关数据情况，包括出生年月、收入、初次来到X市的日期、迁离X市的日期和现在的朋友数量。现假设外来人口的年龄、在X市的居住时间和朋友数量影响他们的收入。试加以证明。

二、要求和代码

一、分析收入的影响因素
#1
#展示数据集的结构
data2 <- read.csv(file="F:/hxpRlanguage/homework2.csv",header=TRUE,sep=",")
str(data2) #显示的结果有一列是多余的，需要删除
data2 <- data2[,1:9]
str(data2) #删完之后的显示效果是正常的没有多余列

#2
#显示前10条数据记录
data2[1:10,]

#3
#将变量名重新命名为英文变量名
cnames <- c("number","birthyear","birthmonth","salary","inyear","inmonth","outyear","outmonth","friends")
colnames(data2) <- cnames
View(data2)

#4
#查找数据集中居住时间小于等于0的异常记录，若存在，从数据集中删除这些异常记录
x2 <- ((data2$outyear-data2$inyear)*12+(data2$outmonth-data2$inmonth))
#View(x2) #①先算出居住时间
data3 <- cbind(data2,x2)
#View(data3) #②使用cbind函数把x2和原数据拼成新的矩阵，方便之后删除异常数据列，并且是127条
list <- which(x2<=0)
data3 <- data3[-list,]
View(data3) #删除异常数据后是125条数据

#5
#展示数据集中因变量与自变量的均值、最小值、中位数、最大值和标准差，要求保留2位小数。
library(lubridate)
date<-Sys.Date() #返回系统当前的时间
nowyear<-year(date) #提取年份
nowmonth<-month(date)  #提取月份
#View(date) #查看现在的日期
#View(month(date)) #查看现在日期中的月份
x1 <- array(1:nrow(data3),dim=c(nrow(data3),1))
for(i in c(1:nrow(data3)) ){
  if(nowmonth-data3[i,"birthmonth"]<0){
   x1[i,1] <- nowyear-data3[i,"birthyear"]-1
  }else{
   x1[i,1] <- nowyear-data3[i,"birthyear"]
  }
}
#View(x1) #算出年龄x1，并加入到数据表中
data4 <- cbind(data3,x1)
View(data4) #加入x1年龄变量的新表展示
x2 <- data4$x2
Mean.x2 <- round(mean(x2),2)
Min.x2 <- round(min(x2),2)
Max.x2 <- round(max(x2),2)
Median.x2 <- round(median(x2),2)
Sd.x2 <- round(sd(x2),2)
cbind(Mean.x2,Min.x2,Max.x2,Median.x2,Sd.x2) #x2居住时间的相关结果
Mean.x1 <- round(mean(x1),2)
Min.x1 <- round(min(x1),2)
Max.x1 <- round(max(x1),2)
Median.x1 <- round(median(x1),2)
Sd.x1 <- round(sd(x1),2)
cbind(Mean.x1,Min.x1,Max.x1,Median.x1,Sd.x1) #x1年龄的相关结果
x3 <- data4$friends
Mean.x3 <- round(mean(x3),2)
Min.x3 <- round(min(x3),2)
Max.x3 <- round(max(x3),2)
Median.x3 <- round(median(x3),2)
Sd.x3 <- round(sd(x3),2)
cbind(Mean.x3,Min.x3,Max.x3,Median.x3,Sd.x3) #x3朋友数量的相关结果
y <- data4$salary
Mean.y <- round(mean(y),2)
Min.y <- round(min(y),2)
Max.y <- round(max(y),2)
Median.y <- round(median(y),2)
Sd.y <- round(sd(y),2)
cbind(Mean.y,Min.y,Max.y,Median.y,Sd.y) #因变量y的相关结果

#6
#计算数据集中因变量和自变量的相关系数，要求保留2位小数。
round(cor(y,x1),2) #y和x1年龄
round(cor(y,x2),2) #y和x2居住时间
round(cor(y,x3),2) #y和x3朋友数量

#7
#分别绘制数据集中因变量与各个自变量的散点图
par(mfrow=c(1,3)) #布局，一行画3个图
plot(x1,y,xlab="年龄x1",ylab="工资y")
plot(x2,y,xlab="居住时间x2",ylab="工资y")
plot(x3,y,xlab="朋友数量x3",ylab="工资y")

#8
#利用多元线性回归模型对数据集中因变量与自变量的关系进行拟合。
lm.xy <- lm(y~x1+x2+x3)
lm.xy
summary(lm.xy) #得到的结果是方程是显著的具有线性关系，但是每一个系数不都是显著的

#9
#对#8中的多元线性回归模型进行诊断，确定异常值记录。
par(mfrow = c(2,2)) #生成四种模型诊断的图形，2行2列
#生成四种模型诊断的图形：①残差与真实值的关系图 ②qq图用来检测其残差是否是正态分布
#③用来检查等方差假设的。在一开始我们的五大假设第二条便是，我们假设预测的模型里方差是一个定值。
#如果方差不是一个定值那么这个模型的可靠性也是大打折扣的。
#④Leverage就是杠杆的意思。这种图的意义在于检查数据分析项目中是否有特别极端的点。
plot(lm)
library(carData)
library(car)
outlierTest(lm.xy) #显示离群点,Bonferroni校正，残差最大的点是136号点

#10
#删除异常值记录后重新利用多元线性回归模型拟合数据。
data4 <- data4[-136,] #删除该点
x1 <- data4$x1
x2 <- data4$x2
x3 <- data4$friends
y <- data4$salary
lm.xy2 <- lm(y~x1+x2+x3) #重新拟合回归模型
lm.xy2

#11
#对#10中的多元线性回归模型进行多重共线性检验，若存在多重共线性，删除相关变量后重新进行拟合。
vif(lm.xy2) #p判断多重共线性0<VIF<10(不存在)

#12
#对#11中的结果进行解释，重点分析年龄、在北京的居住时间和朋友数量如何影响收入。
summary(lm.xy2) #可知年龄和朋友数量对收入有影响，显著性*一颗星

**********************************************************************

二、利用多元线性回归模型预测收入
View(data4) #124条数据
#1
#从数据集中随机抽取50条记录作为预测集，剩下的数据作为训练集。
train0 <- sample(nrow(data4),nrow(data4)-50) #训练集和测试集
trainData <- data4[train0,] #训练数据
testData <- data4[-train0,] #测试数据

#2
#针对训练集，利用多元线性回归模型拟合数据。
lm.xy3 <- lm(trainData[,"salary"]~trainData[,"x1"]+trainData[,"x2"]+trainData[,"friends"])

#3
#对（2）中的多元线性回归模型进行诊断，处理异常值。
summary(lm.xy3)
par(mfrow=c(2,2))
plot(lm.xy3)
outlierTest(lm.xy3)
trainData<-trainData[-c(150,32,82),] #删除异常值，随机的

#4
#对（3）中的多元线性模型进行多重共线性检验并加以处理。
vif(lm.xy3) #p判断多重共线性0<VIF<10(不存在)
salary<-trainData[,"salary"] #引入的数据是训练集的数据
x2<-trainData[,"x2"]
x1<-trainData[,"x1"]
friends<-trainData[,"friends"]
lm.xy3 <-lm(salary~x2+x1+friends)

#5
#针对（4）中的模型，分别利用AIC和BIC选择最优模型。
#AIC检验,赤池信息准则，选择最小的
AIClm<-step(lm.xy3,direction="both")
#BIC检验,贝叶斯信息准则，选择最小的
BIClm<-step(lm.xy3,k=log(nrow(trainData)),direction="both")

#6
#利用预测集进行预测，比较全模型（包含所有自变量）、AIC选择的最优模型、BIC选择的最优模型
#这三个模型预测的准确性大小，并进行解释。
Allmodel<-predict(lm.xy3,testData)
AICmodel<-predict(AIClm,testData)
BICmodel<-predict(BIClm,testData)
#均方误差检验，最小最优，分别计算全模型，AIC,BIC的均方误差
#均方根误差亦称标准误差，均方根误差是预测值与真实值偏差的平方与观测次数n比值的平方根
#标准误差能够很好地反映出测量的精密度
MSE <- function(x){
  mse <- sum((testData[,"salary"]-x)^2)/50
  return(mse)
}
MSE(AICmodel) #AIC/BIC/ALL是误差最小的
MSE(BICmodel)
MSE(Allmodel)

作者: 试试吧 时间: 2021-10-30 17:56
好，谢谢楼主的热情分享的资料

欢迎光临数学建模社区-数学中国 (http://www.madio.net/)