统计学--基于R（第3版）（基于R应用的统计学丛书）作者：贾俊平 习题答案 第八章

8.1

#8.1
load('C:/exercise/ch8/exercise8_1.RData')
exercise8_1
#(1)检验机器对装填量是否有显著影响(α=0.01)
#H0:不显著，H1:显著
#绘制4台机器装填量的箱线图
attach(exercise8_1)
boxplot(填装量~机器,col="gold",main="",ylab="填装量",xlab="机器")
#计算描述统计量
my_summary<-function(x){with(x,data.frame("均值"=mean(填装量),"标准差"=sd(填装量),n=length(填装量)))}
library(plyr)
ddply(exercise8_1,.(机器),my_summary)
#从箱线图可看出，机器3装填量明显低于其他3台，机器1的装填量则相对较高，机器2，4差异不太大，从均值也可看出。
#方差分析表
model_1w<-aov(填装量~机器)
summary(model_1w)
#方差分析模型的参数估计
model_1w$coefficients
#综上所述，Pr(>F)=0.0135，即P>0.01,接受H0,认为机器对装填量没有显著影响
#(2)分析效应量
library(DescTools)
model_1w<-aov(填装量~机器)
EtaSq(model_1w,anova=T)
# eta.sq=0.4070206
#(3)分别采用LSD方法和HSD方法比较哪些机器的装填量之间存在显著差异
#LSD方法
#输出基本结果
library(agricolae)  #3.6.3版本无法下载该package,又下载了4.1.3版本
model_1w<-aov(填装量~机器,data=exercise8_1)
LSD<-LSD.test(model_1w,"机器")
LSD
#由$groups可知，机器1，2，4的平均装填量没有显著差异，机器3的平均装填量与他们三个差异显著
#输出更多信息
library(DescTools)
PostHocTest(model_1w,method="lsd")
#HSD方法
#多重比较的HSD方法
TukeyHSD(model_1w)
#绘制配对差值置信区间的比较图
plot(TukeyHSD(model_1w))
#多重比较的HSD方法
library(agircolae)
HSD<-HSD.test(model_1w,"机器")
HSD
####(3)题分析可照课本221页进行阐述####
#(4)检验装填量是否满足正态性和方差齐性
#正态性检验
#方法1：用Q-Q图检验正态性（直接引用exercise8_1中的数据不行，数据形式不符合要求，需要进行转换）
#绘制每台机器装填量的正态Q-Q图
par(mfrow=c(1,4),cex=0.6,mai=c(0.5,0.5,0.2,0.1))
qqnorm(exercise8_1$填装量[exercise8_1$机器=='机器1'],xlab="期望正态值",ylab="观察值",datax=TRUE,main="机器1的Q-Q图")
qqline(exercise8_1$填装量[exercise8_1$机器=='机器1'],datax=TRUE)
qqnorm(exercise8_1$填装量[exercise8_1$机器=='机器2'],xlab="期望正态值",ylab="观察值",datax=TRUE,main="机器2的Q-Q图")
qqline(exercise8_1$填装量[exercise8_1$机器=='机器2'],datax=TRUE)
qqnorm(exercise8_1$填装量[exercise8_1$机器=='机器3'],xlab="期望正态值",ylab="观察值",datax=TRUE,main="机器3的Q-Q图")
qqline(exercise8_1$填装量[exercise8_1$机器=='机器3'],datax=TRUE)
qqnorm(exercise8_1$填装量[exercise8_1$机器=='机器4'],xlab="期望正态值",ylab="观察值",datax=TRUE,main="机器4的Q-Q图")
qqline(exercise8_1$填装量[exercise8_1$机器=='机器4'],datax=TRUE)
#或将长数据转换为宽数据，再绘制Q-Q图（reshape2中的dcast函数）

#绘制4台机器填装量数据合并后的正态Q-Q图
par(cex=0.8,mai=c(0.7,0.7,0.1,0.1))
qqnorm(exercise8_1$填装量,xlab="期望正态值",ylab="观察值",data=TRUE,main="")
qqline(exercise8_1$填装量,datax=TRUE,col="red",lwd=2)
par(fig=c(0.08,0.5,0.5,0.98),new=TRUE)
hist(exercise8_1$填装量,xlab="填装量",ylab="",freq=FALSE,col="lightblue",cex.axis=0.7,cex.lab=0.7,main="")
lines(density(exercise8_1$填装量),col="red",lwd=2)
box()

#方法2：Shapiro-Wilk和K-S正态性检验
load('C:/exercise8/exercise8_1.RData')
exercise8_1
shapiro.test(exercise8_1$填装量)
ks.test(exercise8_1$填装量,"pnorm",mean(exercise8_1$填装量),sd(exercise8_1$填装量))
#方法1和方法2都可看出，可以认为装填量符合正态分布

#方差齐性检验
#用残差图检验方差齐性
model_1w<-aov(填装量~机器,data=exercise8_1)
par(mfrow=c(1,2),mai=c(0.5,0.5,0.2,0.1),cex=0.6,cex.main=0.7)
plot(model_1w,which=c(1,2))
#Bartlett方差齐性检验
bartlett.test(填装量~机器,data=exercise8_1)
#Levene方差齐性检验
library(car)
leveneTest(填装量~机器,data=exercise8_1)

8.2

#8.2
load("C:/exercise/ch8/exercise8_2.RData")
exercise8_2
#(1)检验管理者的水平不同是否会导致评分的差异显著(α=0.05)，并分析效应量
#H0:不显著，H1:显著
#绘制3层管理者评分的箱线图
attach(exercise8_2)
boxplot(评分~管理者,col="gold",main="",ylab="评分",xlab="管理者")
#箱线图可见，差异显著
#计算描述统计量
my_summary<-function(x){with(x,data.frame("均值"=mean(评分),"标准差"=sd(评分),n=length(评分)))}
library(plyr)
ddply(exercise8_2,.(管理者),my_summary)
#从均值可见，差异显著
#方差分析表
attach(exercise8_2)
model_1w<-aov(评分~管理者)
summary(model_1w)
#Pr(>F)=0.000633，即P<0.05,拒绝H0,认为管理者的水平不同会导致评分的差异显著
#方差分析模型的参数估计
model_1w$coefficients
#效应量分析
library(DescTools)
model_1w<-aov(评分~管理者)
EtaSq(model_1w,anova=T)
#eta.sq=0.6254296
#(2)分别采用LSD方法和HSD方法比较管理者的评分之间的差异
#LSD
#输出基本结果
library(agricolae)
model_1w<-aov(评分~管理者,data=exercise8_2)
LSD<-LSD.test(model_1w,"管理者")
LSD
#由$groups可知，低层管理者与其他二者评分均值差异显著
#输出更多信息
library(DescTools)
PostHocTest(model_1w,method="lsd")
#HSD
#多重比较的HSD方法
TukeyHSD(model_1w)
#绘制配对差值置信区间的比较图
plot(TukeyHSD(model_1w))
#多重比较的HSD方法
library(agricolae)
HSD<-HSD.test(model_1w,"管理者")
HSD
#(3)检验满意度评分是否满足正态性和方差齐性
#正态性检验
#用Q-Q图检验正态性
#绘制每种管理者评分的正态Q-Q图
par(mfrow=c(1,3),cex=0.6,mai=c(0.5,0.5,0.2,0.1))
qqnorm(exercise8_2$评分[exercise8_2$管理者=='高层管理者'],xlab="期望正态值",ylab="观察值",datax=TRUE,main="高层管理者的Q-Q图")
qqline(exercise8_2$评分[exercise8_2$管理者=='高层管理者'],datax=TRUE)
qqnorm(exercise8_2$评分[exercise8_2$管理者=='中层管理者'],xlab="期望正态值",ylab="观察值",datax=TRUE,main="中层管理者的Q-Q图")
qqline(exercise8_2$评分[exercise8_2$管理者=='中层管理者'],datax=TRUE)
qqnorm(exercise8_2$评分[exercise8_2$管理者=='低层管理者'],xlab="期望正态值",ylab="观察值",datax=TRUE,main="低层管理者的Q-Q图")
qqline(exercise8_2$评分[exercise8_2$管理者=='低层管理者'],datax=TRUE)
#绘制3种管理者评分数据合并后的正态Q-Q图
par(cex=0.8,mai=c(0.7,0.7,0.1,0.1))
qqnorm(exercise8_2$评分,xlab="期望正态值",ylab="观察值",data=TRUE,main="")
qqline(exercise8_2$评分,datax=TRUE,col="red",lwd=2)
par(fig=c(0.08,0.5,0.5,0.98),new=TRUE)
hist(exercise8_2$评分,xlab="评分",ylab="",freq=FALSE,col="lightblue",cex.axis=0.7,cex.lab=0.7,main="")
lines(density(exercise8_2$评分),col="red",lwd=2)
box()
#Shapiro-Wilk正态性检验
shapiro.test(exercise8_2$评分)
#K-S正态性检验
ks.test(exercise8_2$评分,"pnorm",mean(exercise8_2$评分),sd(exercise8_2$评分))
#方差齐性检验
#用残差图检验方差齐性
model_1w<-aov(评分~管理者,data=exercise8_2)
par(mfrow=c(1,2),mai=c(0.5,0.5,0.2,0.1),cex=0.6,cex.main=0.7)
plot(model_1w,which=c(1,2))
#Bartlett方差齐性检验
bartlett.test(评分~管理者,data=exercise8_2)
#Levene方差齐性检验
library(car)
leveneTest(评分~管理者,data=exercise8_2)
#(4)假定不同层次管理者的满意度评分不满足方差分析的假定，采用Kruskal-Wallis检验进行分析
attach(exercise8_2)
kruskal.test(评分~管理者)

8.3

#8.3
load("C:/exercise/ch8/exercise8_3.RData")
exercise8_3
#(1)分析3个企业生产的电池平均寿命之间有无显著差异(α=0.05)，并分析效应量
#H0:不显著，H1:显著
#绘制3个企业电池平均寿命的箱线图
attach(exercise8_3)
boxplot(电池寿命~企业,col="gold",main="",ylab="电池寿命",xlab="企业")
#箱线图可见B与其他两个差异显著
#计算描述统计量
my_summary<-function(x){with(x,data.frame("均值"=mean(电池寿命),"标准差"=sd(电池寿命),n=length(电池寿命)))}
library(plyr)
ddply(exercise8_3,.(企业),my_summary)
#均值也可见B与其他两个差异显著
#方差分析表
attach(exercise8_3)
model_1w<-aov(电池寿命~企业)
summary(model_1w)
#Pr(>F)=0.00031，即P<0.05,拒绝H0,
#方差分析模型的参数估计
model_1w$coefficients
#绘制均值图
library(gplots)
plotmeans(电池寿命~企业,data=exercise8_3)
#(2)如果有差异，用HSD方法检验哪些企业之间有差异
#多重比较的HSD方法
TukeyHSD(model_1w)
#绘制配对差值置信区间的比较图
plot(TukeyHSD(model_1w))
#多重比较的HSD方法
library(agricolae)
HSD<-HSD.test(model_1w,"企业")
HSD
#B与A,C有差异

8.4

#8.4
load("C:/exercise/ch8/exercise8_4.RData")
exercise8_4
#(1)检验不同品种和施肥方案对产量的影响是否显著(α=0.05)
#绘制品种和施肥方案的箱线图
attach(exercise8_4)
boxplot(收获量~品种+施肥方案,col=c("gold","green","red"),ylab="收获量",xlab="品种与施肥方案")
#按品种和施肥方案交叉分类计算均值和标准差
library(reshape)
library(agricolae)
mystats<-function(x)(c(mean=mean(x),sd=sd(x)))#将函数中n=length(x)去掉，因为收获量是浮点数，不能计算其有效长度。
dfm<-melt(exercise8_4,measure.vars="收获量",id.vars=c("品种","施肥方案"))
cast(dfm,品种+施肥方案+variable~.,mystats)
#主效应方差分析
model_2wm<-aov(收获量~品种+施肥方案)
summary(model_2wm)
#由于p1=0.00332,p2=0.00195,均小于0.05，所以认为不同品种和施肥方案对产量的影响均显著
#主效应方差分析模型的参数估计
model_2wm$coefficients
#(2)计算偏效应量并进行分析
library(DescTools)
EtaSq(model_2wm,anova=T)
#eta.sq.part列为偏效应量，结果表明，在排除施肥方案的影响后，品种因子解释了产量误差的70.70231%；再排除品种的影响后，施肥方案解释了产量误差的69.70766%

8.5

#8.5
load("C:/exercise/ch8/exercise8_5.RData")
exercise8_5
#(1)分析路段、时段以及路段和时段的交互作用对行车时间的影响(α=0.05)
#交互效应方差分析表
model_2wi<-aov(行车时间~时段+路段+时段:路段,data=exercise8_5)
summary(model_2wi)
#结果显示，时段和路段的p均小于0.05，表示两个因子对行车时间的影响显著，而检验交互效应的p=0.997>0.05,表示交互效应对行车时间的影响不显著
#交互效应方差分析模型的参数估计
model_2wi$coefficients
#绘制时段和路段的主效应和交互效应图
library(HH)
interaction2wt(行车时间~时段+路段,data=exercise8_5)
#(2)计算偏效应量并进行分析
model_2wi<-aov(行车时间~时段+路段+时段:路段,data=exercise8_5)
library(DescTools)
EtaSq(model_2wi,anova=T)
#eta.sq.part列为偏效应量，结果表明，在排除路段的影响后，时段因子解释了行车时间的74.55777387%；在排除时段的影响后，路段因子解释了行车时间的64.72503358%

8.6

#8.6
load("C:/exercise/ch8/exercise8_6.RData")
exercise8_6
#(1)检验广告方案、广告媒体及其交互作用对销售量的影响是否显著(α=0.05)
#交互效应方差分析表
model_2wi<-aov(销售量~广告方案+广告媒体+广告方案:广告媒体,data=exercise8_6)
summary(model_2wi)
#结果显示，广告方案的p<0.05,表明广告方案对销售量的影响显著，而广告媒体以及交互作用的p均大于0.05，表明它们对销售量的影响均不显著
#交互效应方差分析模型的参数估计
model_2wi$coefficients
#绘制广告方案与广告媒体的主效应和交互效应图
library(HH)
interaction2wt(销售量~广告方案+广告媒体,data=exercise8_6)
#(2)计算偏效应量并进行分析
model_2wi<-aov(销售量~广告方案+广告媒体+广告方案:广告媒体,data=exercise8_6)
library(DescTools)
EtaSq(model_2wi,anova=T)
#eta.sq.part列为偏效应量，结果表明，在排除广告媒体的影响后，广告方案解释了销售量的78.18182%；在排除广告方案的影响后，广告媒体解释了销售量的33.33333%