多元统计分析及R语言建模(第五版)部分课后习题代码演示

前言:我是lst,这是本学期课程内容,根据上课内容、参考书本、CSDN社区等来源,完成的部分课程习题代码。如有问题,欢迎大家批评指正


library(openxlsx)                                     #加载读取Excel数据包

#【输出设置】
#setwd("C:/Users/lst89/Documents/mvexer5") #设置目录
options(digits=4)
par(mar=c(4,4,2,1))

#第二章p57-2-1
R=matrix(c(1,0.8,0.26,0.67,0.34,0.8,1,0.33,0.59,0.34,0.26,0.33,1,0.37,0.21,0.67,0.59,0.37,1,0.35,0.34,0.34,0.21,0.35,1),nrow = 5,ncol = 5);
R                                                              #输入数据
solve(R)                                                   #求逆矩阵
R.e=eigen(R,symmetric=T)                     #symmetric是判断是否为对称阵,
R.e                                                           #求矩阵的特诊值
R.e $ vectors%*%diag(R.e $ values)%*%t(R.e $ vectors)#特征向量

#第二章p57-2-2
library(openxlsx)                                      #加载读取Excel数据包
E2.2=read.xlsx('mvexer5.xlsx','E2.2');    
E2.2                                                         #读取mvexer5.xlsx表格E2.2数据
breaks = seq(0,3000,by = 300)               #按组距为300编制频数表
breaks
hist(E2.2 $ X,breaks,col = 1:7,xlab = "工资(元)",ylab = "频数")#以工资x为横轴,频数y为纵轴,将数据划分为0-3000并以300为度量,绘制7列的彩色直方图
hist(E2.2 $ X ,breaks,freq = F,col = 1:7,xlab = "工资(元)",ylab = "频率")
Cumsum <- cumsum(E2.2 $ X)
cumsum
M <- seq(0,96000,by = 3000)

hist(Cumsum,M,freq = F,col = 1:12,las = 3,xlab = "工资(元)",ylab = "累积频率")#绘制出累计频率直方图
H = hist(E2.2 $ X,breaks = seq(900,3000,300))#正态概率图
names(H)
data.frame('组中距' = H $ mids,'频数' = H $ counts,'频率' = H $ density*300,'累积频率' = cumsum(H $ density*300))#


#第二章p57-2-3
library(openxlsx)                                             #加载读取Excel数据包
E2.3=read.xlsx('mvexer5.xlsx','E2.3');          
E2.3                                                                #读取mvexer5.xlsx表格E2.2数据
str(E2.3)                                       
summary(E2.3)                                                #对数据进行基本统计分析


#第三章P84-2.1
library(openxlsx)
E3.2 = read.xlsx('mvexer5.xlsx',sheet = 'E3.2',rowNames = TRUE)          
#设定参数rowNames=TRUE,即可将第一列字符变量变成数据框的行名,供后期使用
E3.2
#在Excel文件中mvexer5.xlsx的表单d3.2中选择A1:E22,并复制到剪切板
dat = read.table("clipboard",header = T)               #将剪切板数据读入数据框dat中
dat
#数据框标记转换函数
msa.X <- function(df){                                          #将数据框第一列设置为数据框行名
  X = df[,-1]                                                          #删除数据框df的第一列并赋给X
  rownames(X) = df[,1]                                       #将df的第一列值赋给X的行名
  X                                                                      #返回新的数值数据框=return(X)
}
E3.2 = msa.X(dat)
E3.2
barplot(apply(E3.2,2,mean))                              #按行作均值条形图
barplot(apply(E3.2,1,mean),las = 3)                  #修改横坐标标记
barplot(apply(E3.2,2,mean))                             #按列作均值条图
barplot(apply(E3.2,2,median))                          #按列作中位数条图
barplot(apply(E3.2,2,median),col = 1:8)           #按列取色
boxplot(E3.2)                                                    #按列作箱尾图
boxplot(E3.2,horizontal = T)                             #箱尾图中图形按水平放置


#四p119-2-1
library(openxlsx)                                               #加载读取Excel数据包
E4.1=read.table("clipboard",header = T)
E4.1

plot(x,y,main = '散点图',xlab = '每周加班时间(小时)',ylab = '每周签发的新保单数目(张)')                                                             #绘制散点图

cor(E4.1)                                                          #相关系数
lm4.1 <- lm(E4.1)
lm4.1

#估计值
square_sigma <- t(E4.1)/(10-1-1)#square_sigma <- t(x_hat - y)%*%(x_hat - y)/(10-1-1)
square_sigma                                    

y = c(3.5,1,4,2,1,3,4.5,1.5,3,5)
x = c(825,215,1070,550,480,920,1350,325,670,1215)
y_hat <- 46.15 + 251.17*y
s <- t(y_hat - x)%*%(y_hat - x)/(10-1-1)
s
(summary(lm4.1) $ s)^2

#求方差分析
SR <- t(y_hat - mean(x))%*%(y_hat - mean(x))
ST <- t(x - mean(x))%*%(x - mean(x))
s_R <- SR/ST
s_R                                       
(summary(lm4.1) $ r.squared)
anova(lm4.1)                              

#对回归方程作残差图分析
res <- residuals(lm4.1)
res
plot(y,res,main='残差散点图',xlab='每周签发的新保单数目',ylab='残差')
plot(lm4.1)            

#计算1000张要加班的时间
lm4.1_1 <- lm(x ~ y,data = ee4.1)
predict(lm4.1_1,newdata = data.frame(y = 1000))
lm4.1_1 <- lm(y ~ x,data = ee4.1)
predict(lm4.1_1,newdata = data.frame(x = 1000))     


#四p119-2-2
library(openxlsx)
E4.2 = read.xlsx('mvexer5.xlsx',sheet = 'E4.2',rowNames = T)
(lm4.2 = lm(y ~ x1 + x2,data = E4.2))                 #显示多元线性回归模型


library(openxlsx)                                                 #加载读取Excel数据包
E4.2=read.table("clipboard",header = T)
E4.2
(lm4.2=lm(y~x1+x2,data=E4.2))
coef.sd <- function(lm4.2){                                 #标准回归系数
  b = lm4.2 $ coef
  b
  si = apply(lm4.2 $ model,2,sd)
  si
  bs = b[-1] * si[-1]/si[1]
  bs
}

coef.sd(lm4.2)
                                                  #在5%的显著水平下确定每一解释变量与依赖变量是否呈线性关系
summary(lm4.2)
                                                  #计算简单相关系数和复相关系数
cor(E4.2)


#四119-2-3
#回归模型并解释各系数
E4.3 = read.table("clipboard",header = T)
E4.3
#确定学生的GPA和年龄是否能真正用来解释起始工资样本的变化
(lm4.3 = lm(起始工资~ GPA + 年龄,data = E4.3))
summary(lm4.3)
R4.3 = summary(lm4.3) $ r.sq
R4.3
R_2_4.3 = sqrt(R4.3)
R_2_4.3
#预测某GPA为3.00,年龄为24岁的毕业生的起始工资
predict(lm4.3,newdata = data.frame(GPA = 3,年龄 = 24))

#四120-2-4
#计算y,x1,x2,x3的相关系数矩阵并绘制矩阵散点图

library(openxlsx)
E4.4 = read.xlsx('mvexer5.xlsx',sheet = 'E4.4')
cor4.4 <- cor(E4.4[,-1])      #去除第一列编号数据后剩余的样本计算相关系数矩阵
cor4.4

pairs(E4.4)#绘制散点图

#求y关于x1,x2,x3的多元线性回归方程
lm4.4 <- lm(y ~ x1 + x2 + x3,data = E4.4)      #建立回归方程
summary(lm4.4)
#所求得的方程做拟合优度检验
summary(lm4.4)
(R4.4 = summary(lm4.4)$r.sq)
(R_2_4.4 = sqrt(R4.4))

#对回归方程做显著性检验,对每一个回归系数做显著性检验
summary(lm4.4)
#如果有的回归系数没通过显著性检验,将其剔除,重新建立回归方程,再做回归方程的显著性检验和回归系数的显著性检验
lm4.4_drop3 <- lm(y ~ x1 + x2,data = E4.4)
lm4.4_drop3
summary(lm4.4_drop3)

#对每一个回归系数做显著性检验
summary(lm4.4_drop3)
t1 <- 1.942
t2 <- 2.465
t3 <- 1.178
qt(1-0.05,7)

#使用逐步回归分析的逐步筛选方法获得一个最优的回归模型
lm.step = step(lm4.4,direction = 'forward')           #向前引入法变量选择结果

lm.step = step(lm4.4,direction = 'backward')       #向后剔除法变量选择结果

lm.step = step(lm4.4,direction = 'both')               #逐步筛选法变量选择结果


#五p144-2-1
library(openxlsx)
E5.1 = read.xlsx('mvexer5.xlsx',sheet = 'E5.1')
E5.1
anova(lm(S ~ factor(F),data = E5.1))

#五p145-2-2
library(openxlsx)
E5.2 = read.xlsx('mvexer5.xlsx',sheet = 'E5.2')
E5.2
summary(E5.2)
anova(lm(P ~ factor(A)+factor(B),data=E5.2))

#五p145-2-3
library(openxlsx)
E5.3 <- read.xlsx('mvexer5.xlsx',sheet = 'E5.3')
E5.3
summary(E5.3)
logit.glm = glm(G ~ x1 + x2,data = E5.3)
summary(logit.glm)

logit.step <- step(logit.glm,direction = 'both')    #逐步筛选法变量选择
summary(logit.step)
pre1 <- predict(logit.glm,data.frame(x1 = 131,x2 = -2))
pre1

#第七章p212-2
library(openxlsx)
E7.2=read.xlsx("mvexer5.xlsx",sheet = 'E7.2',rowNames = T)
E7.2
plot(E7.2,gap=0)                                   #绘制散点图
D =dist(E7.2)                                        #计算距离矩阵
D
hc = hclust(D,"single")                          #聚类 最短距离法
hc
names(hc)
data.frame(hc $ merge,hc $ height)
plot(hc)

hc = hclust(D)                                      #聚类 最长距离法
hc
names(hc)
data.frame(hc $ merge,hc $ height)
plot(hc)

plot(hclust(D,"median"))                     #聚类 中间距离法
plot(hclust(D,"average"))
plot(hclust(D,"centroid"))
plot(hclust(D,"ward.D"))
plot(hclust(D,"ward.D2"))
plot(hc)
rect.hclust(hc,k = 4,border = 'pink')    #综上 认为平均距离好。分四类

  • 12
    点赞
  • 128
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值