【高级数理统计R语言学习】5 协方差分析

一、背景
数据集展示了已迁离X市的高学历外来人口现在的月收入、性别、迁入X市和迁离X市的日期、教育程度和职业这些方面的数据。试分析性别、在X市的居住时间、教育程度和职业对外来人口的收入是否有显著影响以及有怎样的影响。要求分析教育程度与在X市的居住时间、职业与在X市的居住时间的交互作用对收入的影响。

二、要求和代码

#**************************变量关系问题******************************
#1
#利用R读取数据，将变量名重新命名。
data5 <- read.csv(file="F:/hxpRlanguage/homework5.csv",header=TRUE,sep=",",stringsAsFactors = F)
data5 <- na.omit(data5) #删除空行
cnames <- c("number","Gender","inyear","inmonth","outyear","outmonth","Education","Career","Income")
colnames(data5) <- cnames#2
#显示数据集中变量的属性情况。
str(data5)#3
#显示数据集的前10条记录。
data5[1:10,]#4
#利用R编写程序对性别、教育程度和职业变量进行归类，不能用Excel处理数据。
#类别划分标准见Word文件“数据编码要求”。性别的变量名为Gender，取值为0、1；教育程度的变量名为Education，取值为0、1
#职业的变量名为Career，取值为1、2、3、4、5。注意：样本中有填写的半结构化数据。
#①对性别进行归类
data5$fGender[data5$Gender=="A. 男"]<-"1"
data5$fGender[data5$Gender=="B. 女"]<-"0"
#②对教育程度进行归类
data5$fEducation[data5$Education=="E. 大专" | data5$Education=="F. 大学本科"]<-"0"
data5$fEducation[data5$Education=="G. 硕士" | data5$Education=="H. 博士"]<-"1"
#③对职业进行归类
data5$fCareer[substring(data5$Career,1,2)=='O1'|data5$Career=="企业/公司厂长经理或管理者"|data5$Career=="党政机关领导人"|data5$Career=="私人企业主"|data5$Career=="企业党委书记"|data5$Career=="外商代理人"|data5$Career=="房产商"]<-1
data5$fCareer[substring(data5$Career,1,2)=='O2'|data5$Career=="工商税务干部"|data5$Career=="公安政法干部"|data5$Career=="银行职员"|data5$Career=="外资公司职员"|data5$Career=="外贸公司职员"|data5$Career=="机关事业单位一般职员"|data5$Career=="民政工作者"|data5$Career=="秘书"]<-2
data5$fCareer[substring(data5$Career,1,2)=='O3'|data5$Career=="律师"|data5$Career=="节目主持人"|data5$Career=="歌星"|data5$Career=="医生"|data5$Career=="大学教师"|data5$Career=="科学家"|data5$Career=="演员"|data5$Career=="音乐家"|data5$Career=="画家"|data5$Career=="记者"|data5$Career=="科研人员"|data5$Career=="工程师"|data5$Career=="会计"|data5$Career=="中小学教师"]<-3
data5$fCareer[substring(data5$Career,1,2)=='O4'|data5$Career=="经纪人"|data5$Career=="企业/公司一般职员"|data5$Career=="物资管理干部"|data5$Career=="个体户/自由职业者"|data5$Career=="销售员"|data5$Career=="营业员"]<-4
data5$fCareer[substring(data5$Career,1,2)=='O5'|substring(data5$Career,1,2)=='O6'|substring(data5$Career,1,2)=='O7'|data5$Career=="远洋轮船员"|data5$Career=="出租车司机"|data5$Career=="公交司机"|data5$Career=="火车司机"|data5$Career=="电工"|data5$Career=="建筑工人"|data5$Career=="纺织工人"|data5$Career=="家电维修工"|data5$Career=="厨师"|data5$Career=="邮递员"|data5$Career=="快递员"|data5$Career=="宾馆饭店服务员"|data5$Career=="理发员"|data5$Career=="保育员"|data5$Career=="保姆"|data5$Career=="A清洁工人"|data5$Career=="勤杂工"|data5$Career=="农民"]<-5#5
#性别、教育程度和职业变量归类完毕后，将它们转换成因子变量类型。
data5$fGender <- as.factor(data5$fGender)
data5$fEducation <- as.factor(data5$fEducation)
data5$fCareer <- as.factor(data5$fCareer)
str(data5)#6
#利用R编写程序计算个体在北京的居住时间，单位为月，变量名为Duration。不能用Excel处理数据。
Duration <- (data5$outyear-data5$inyear)*12+(data5$outmonth-data5$inmonth)
data5 <- cbind(data5,Duration)#7
#删除居住时间小于6个月的数据记录。
dl <- which(Duration<6)
data5 <- data5[-dl,]#8
#显示数据框中变量的属性。
str(data5)#9
#分别统计性别、教育程度和职业变量每个水平的个数。
table(data5$fGender)
table(data5$fEducation)
table(data5$fCareer)#10
#职业类型5的观察值太少，将其从总样本中删除。
data5 <- data5[-c(which(data5[,"fCareer"]=="5")),]
data5$logIncome<-log(data5$Income)  #转对数收入
data5.11 <- data5#11
#显示更新后的数据集的前10条记录。
data5[1:10,]#12
#分别画图展示数据集中对数收入对性别、教育程度、职业和居住时间的关系。
par(mfrow = c(2, 2))
boxplot(data5$logIncome~data5$fGender,xlab="性别",ylab="对数收入")
boxplot(data5$logIncome~data5$fEducation,xlab="教育水平",ylab="对数收入")
boxplot(data5$logIncome~data5$fCareer,xlab="职业",ylab="对数收入")
plot(data5[,"logIncome"]~data5[,"Duration"],xlab="居住时间",ylab="对数收入")#13
#作对数收入对性别、教育程度、职业、居住时间的协方差分析，检验职业与居住时间的交互作用是否存在。
#lm5.0 <- lm(logIncome~fGender+fEducation+fCareer+Duration+fCareer*Duration,data=data5)
#anova(lm5.0)
#summary(lm5.0)
lm5.1 <- lm(logIncome~fGender+fEducation+fCareer+Duration,data=data5)
anova(lm5.1)
summary(lm5.1)  #消失的参照物是gender0+edu0+career1=女性低教育程度领导干部为主的群体
lm5.2 <- lm(logIncome~fCareer*Duration,data=data5)  #消失的参照物是career1领导干部
anova(lm5.2)
summary(lm5.2)  #Career4是以商业为主的群体#14
#对数据集中的变量Duration进行归类，类别的划分标准见Word文件“数据编码要求”，归类后的变量取名为fDuration，取值为1、2、3。
#View(data5)
fDuration<-array(1:nrow(data5),dim= c(nrow(data5),1))
for (i in c(1:nrow(data5)) ) {if (data5[i,"Duration"]>=6&data5[i,"Duration"]<=24){fDuration[i,1]<- "1"}else if(data5[i,"Duration"]<=84&data5[i,"Duration"]>24){                              fDuration[i,1]<- "2"}else{fDuration[i,1]<- "3"}
}
names(fDuration)<-"fDuration"
data5 <- cbind(data5,fDuration)
str(data5)#15
#统计变量fDuration各个水平的个数。
table(data5$fDuration)#16
#画图展示数据集中对数收入与fDuration的关系。
par(mfrow = c(1, 1))
boxplot(data5$logIncome~data5$fDuration,xlab="居住时间",ylab="对数收入")#17
#分别对性别、教育程度、职业、fDuration各因素进行方差齐性检验。
library("carData")
library("car")
leveneTest(logIncome~fGender,data=data5)
leveneTest(logIncome~fEducation,data=data5)
leveneTest(logIncome~fCareer,data=data5)
leveneTest(logIncome~fDuration,data=data5)#18
#作对数收入对性别、教育程度、职业、fDuration的多因素方差分析，检验职业与fDuration的交互作用是否存在。
lm5.3 <- lm(logIncome~fGender+fEducation+fCareer+fDuration,data=data5)
anova(lm5.3)
summary(lm5.3) #消失的参照组是gender0+edu0+career1+dura1=女性教育程度低的以领导干部为主的群体6-24月
#检验职业与fDuration的交互作用
lm5.4 <- lm(logIncome~fCareer*fDuration,data=data5)
anova(lm5.4)
summary(lm5.4) #消失的是career1+dura1=领导干部居住6-24个月#19
#将（18）中的方差分析模型与（13）中的协方差分析模型进行对比，分析两者结果的差异。
anova(lm5.3) #方差模型
anova(lm5.1) #斜方差模型#***********************预测问题*******************************************
#1
#将（11）中数据集分为模型检验集和预测集两部分，随机抽取50条记录作为预测集。
sub1 <- sample(nrow(data5.11),50,replace=F)  #共212条数据，随机抽取50条记录作为预测集
data5.11m <- data5.11[-sub1,]  #模型检验集
data5.11p <- data5.11[sub1,]   #预测集#2
#对检验集进行协方差分析并对模型进行诊断。
#检验变量是否存在显著性差异；summary来分析各个变量对对数收入的影响
lm5.5 <- lm(logIncome~fGender+fEducation+fCareer+Duration,data=data5.11m)
anova(lm5.5) #方差分析表
summary(lm5.5) #参数估计表,参照组gender0+edu0+career1=女性低教育领导干部
par(mfrow=c(2,2)) #模型诊断,qq图的点大致分布在直线两侧，残差服从正态分布的假设成立
plot(lm5.5,which=c(1:4))#3
#删除异常值后重新进行协方差分析。
data5.11mr <- data5.11m[-158,] #删除异常值
lm5.6 <- lm(logIncome~fGender+fEducation+fCareer+Duration,data=data5.11mr)
anova(lm5.6,test="Chisq") #方差分析表
summary(lm5.6) #参数估计表#4
#分别以AIC和BIC标准选择最优预测模型。
lm.aic=step(lm5.6,trace=F)
anova(lm.aic,test="Chisq")
#summary(lm.aic)
lm.bic=step(lm5.6,k=log(nrow(data5.11mr)),trace=F)
anova(lm.bic,test="Chisq")
#summary(lm.bic)#5
#利用最优预测模型对预测集进行预测，并计算出预测精度。
#使用平均相对误差的方法
logIncome.hat=predict(lm.aic,data5.11p)
r0=exp(data5.11p$logIncome) #真实值收入
#kan <- cbind(logIncome.hat,data5.11p$logIncome)
#View(kan)
#kan1 <- exp(logIncome.hat)
#kan2 <- exp(data5.11p$logIncome)
#kankan <- cbind(kan1,kan2)
#View(kankan)
#kan3 <- abs(kan1-kan2)
#View(kan3)
#kan4 <- kan3/kan2
#View(kan4)
r1=exp(logIncome.hat)-exp(data5.11p$logIncome) #预测值-真实值
r2=abs(r1) #取绝对值
r3=r2/r0 #差的绝对值/真实值
mean(r2)
mean(r3)#添加一个预测值的对比，AIC和BIC，对数收入
y.aic=predict(lm.aic,data5.11p) #aic预测值
y.bic=predict(lm.bic,data5.11p) #bic预测值
y0=data5.11p$logIncome #真实值
y1=(y0-y.aic)/y0
y2=(y0-y.bic)/y0
perror <- abs(as.data.frame(cbind(y1,y2)))
sapply(perror,mean)

【高级数理统计R语言学习】5 协方差分析相关推荐

【高级数理统计R语言学习】9 无序多分类分析
一.背景数据集展示了人们休闲的相关数据.试分析年龄.性别.教育程度.月收入对人们的休闲方式是否有显著影响以及有怎样的影响. 二.要求和代码 #1 #利用R读取数据 data9 <- read. ...
【高级数理统计R语言学习】7 定序回归
一.背景数据集展示了用户使用微博的基本情况,包括参与微博社区的层次,用户的年龄.性别.教育程度.月收入和使用微博的时间,试分析这些变量对用户参与微博社区的层次有什么样的影响?同时,对用户参与微博社区 ...
R语言学习笔记——高级篇：第十四章-主成分分析和因子分析
R语言 R语言学习笔记--高级篇:第十四章-主成分分析和因子分析文章目录 R语言前言一.R中的主成分和因子分析二.主成分分析 2.1.判断主成分的个数 2.2.提取主成分 2.3.主成分旋转 ...
R语言学习手记 (1)
R语言学习手记 (1) 经管的会计和财管都会学数据统计与分析R语言这门课,加上我也有点兴趣,就提前选了这门课,以下的笔记由老师上课的PPT.<R语言编程艺术>和<R语言数据科学> ...
R语言学习笔记（1~3）
R语言学习笔记(1~3) 一.R语言介绍 x <- rnorm(5) 创建了一个名为x的向量对象,它包含5个来自标准正态分布的随机偏差. 1.1 注释由符号#开头. #函数c()以向量的形式输 ...
R语言学习之数据分析实战(一)
R语言学习之数据分析实战(一) 一.线性回归回归(regression):通常指那些用一个或多个预测变量,也称自变量或解释变量,来预测响应变量,也称为因变量.效标变量或结果变量的方法. 普通最小二乘 ...
R语言学习笔记——入门篇：第三章-图形初阶
R语言 R语言学习笔记--入门篇:第三章-图形初阶文章目录 R语言一.使用图形 1.1.基础绘图函数:plot( ) 1.2.图形控制函数:dev( ) 补充--直方图函数:hist( ) 补充- ...
R语言回归模型协方差分析（Analysis of Covariance）
R语言回归模型协方差分析(Analysis of Covariance) 目录 R语言回归模型协方差分析(Analysis of Covariance) 输入数据 ANCOVA分析
r语言c函数怎么用,R语言学习笔记——C#中如何使用R语言setwd()函数
在R语言编译器中,设置当前工作文件夹可以用setwd()函数. > setwd("e://桌面//") > setwd("e:\桌面\") > ...

【高级数理统计R语言学习】5 协方差分析

【高级数理统计R语言学习】5 协方差分析相关推荐

最新文章

热门文章